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ABSTRACT 

A  laminated  shallow  shell  approach  that  includes  von  Karman  geometric 
nonlinearity  and  parabolic  transverse  shear  deformation  is  posed  in  differential  operator 
form.  Trigonometric  series  are  assumed  for  each  of  the  five  shell  displacement  degrees  of 
fireedom  for  the  subsequent  nonlinear  Galerkin  solution  resulting  in  5n^  shnultaneous 
algebraic  equations  where  n  is  the  number  of  displacement  terms  assumed  in  the  series. 
The  Galerkin  nonlinear  solution  is  computationally  intensive.  The  response  of  several 
laminate  geometries  subjected  to  transverse  loadings  are  found.  Thicker  plates  and  shells 
generally  exhibit  more  flexible  nondimensional  displacement  response  compared  to  thinner 
geometries  in  both  linear  and  nonlinear  analyses  due  to  transverse  shear  deformation.  The 
nondimensional  shell  response  is  examined  by  using  the  Batdorf-Stein  shell  parameter  for 
laminated  constructions.  Quasi-isotropic  shallow  shells  undergo  significant  transverse 
shear  flexibility  in  the  thicker  geometries  as  given  by  the  nondimensional  shell  crown 
deflection.  However,  the  nondimensional  crown  deflection  in  the  deeper  shell  response  is 
not  much  influenced  by  sheU  thickness.  For  flat  plates,  geometric  nonlinearity  lessens  the 
influence  of  transverse  shear  flexibility  when  compared  to  linear  solutions  due  to 
membrane  stretching  resistance. 
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I.  Introduction 


The  increasingly  higher  performance  demanded  of  advanced  composite  materials, 
especially  in  weight  optimized  structures,  requires  more  accurate  means  of  analysis  and 
design.  Traditional  techniques  for  thin-waUed  aerospace  structures  include  classical 
Kirchhoff  plate  and  Love  shell  theories.  These  theories  neglect  the  transverse  stress 
components  and  therefore  are  inadequate  in  many  cases  due  to  the  sensitivity  in  laminated 
materials  to  transverse  shear  deformation  (Whitney  and  Pagano,  1970).  Consequently, 
higher  order  theories  have  been  developed  that  relax  the  requirement  in  the  classical 
approaches  where  normals  to  the  shell  middle  plane  before  deformation  remain  normal 
throughout  deformation;  see  for  example  Lo  et  al  (1977),  Kwon  and  Akin  (1987), 

Whitney  and  Sun  (1974),  and  the  review  articles  Kapania  (1989),  Kapania  and  Reciti 
(1989),  and  Noor  and  Burton  (1989, 1990).  Generally,  Reissner-Mindlin  (R-M) 

(Reissner,  1945,  Mindlin,  1951)  shear  deformation  theory  requires  plane  sections  remain 
plane,  giving  a  constant  transverse  shear  through  the  shell  thickness.  R-M  plate  and  shell 
approaches  then  require  shear  correction  factors  that  depend  on  the  specific  laminate 
(Whitney,  1973)  or  the  boundary  conditions  and  loading,  as  described  by  Pai  (1995). 

More  recently,  a  parabolic  distribution  of  transverse  shear  through  the  thickness  has  been 
proposed  by  Levinson  (1980),  Murthy  (1981),  and  Reddy  (1984)  for  plates  and  Reddy  and 
Liu  (1985)  for  shells.  This  approach  gives  five  coupled  nonlinear  differential  equations  as 
in  the  R-M  approaches.  Yet,  for  only  a  small  additional  complexity,  added  accuracy  is 
gained  without  the  ambiguities  of  shear  correction  factors.  Exact  solutions  to  any  of  these 
shell  approaches  are  rare;  consequently,  approximate  techniques  are  widely  used. 

Numerical  solutions  to  Reissner-Mindlin  theory  based  on  finite  element  methods 
are  susceptible  to  shear  locking  of  thin  plate/shell  meshes  unless  reduced  integration  of  the 
stiffness  terms  (Zienkiewicz  et  al,  1971,  Parisch,  1979 )  or  other  numerical  fix  is 
performed.  (Kui,  1985,  Park  and  Stanley,  1986)  The  higher  order  parabolic  shear  theory 
including  geometric  nonlinearity  has  been  solved  via  finite  element  techniques  for 
laminated  plates  by  Putcha  and  Reddy  (1986)  and  laminated  shells  by  Tsai  et  al  (1991). 
Finite  element  meshes  based  on  the  parabolic  shear  approach  do  not  shear  lock  for  a  wide 
range  of  plate/shell  geometries. 

A  related  solution  approach  to  the  finite  element  method  is  the  modified  Galerkin 
technique.  Giri  and  Simitses  (1980)  used  this  technique  to  solve  geometrically  nonlinear 
plate  equations  based  on  classical  assumptions.  \Vhitney  (1984)  gave  buckling  loads  for 
shallow  laminated  panels  also  neglecting  transverse  shear  deformation.  Xu  and  Shen 
(1986)  used  the  Galerkin  method  to  approximate  the  classical  nonlinear  response  of 
laminated  plates  by  assuming  quintic  B-spline  functions.  Bowlus  et  al  (1987)  used  the 
Galerkin  technique  to  examine  the  effect  of  transverse  shear  flexibility  and  rotary  inertia 
assuming  the  R-M  approach.  Palazotto  and  Linneman  (1991)  found  buckling  loads  for 
laminated  shells  assuming  Sanders  equations  and  parabolic  transverse  shear  distributions 
through  the  shell  thickness.  The  influence  of  transverse  shear  and  radius  of  curvature  is 
examined,  Tighe  and  Palazotto  (1994)  assumed  Sanders  shell  equations  with  nonzero 
transverse  normal  strain. 

The  present  study  solves  laminated  shallow  shell  equations  including  parabolic 
transverse  shear  deformation  in  the  presence  of  geometric  nonlinearity  by  the  modified 
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Galerkin  technique  (Sokolnikoff,  1956).  The  approach  is  &st  cast  into  operator  form, 
resulting  in  five  coupled  nonlinear  differential  equations.  Next,  a  double  Fourier  series  is 
assumed  for  each  of  the  displacements  and  the  Galerkin  equations  are  found  for  simple 
support  boundary  conditions.  Finally,  the  response  of  numerous  plate  and  shell  structures, 
both  linear  and  nonlinear,  are  presented. 

II.  Theory 

The  present  study  will  approximate  the  response  of  a  shallow  cylindrical  shell 
geometry  under  transverse  loading,  see  Figure  1.  Assume  the  following: 

1.  the  shell  is  thin  such  that  an  approximate  state  of  plane  stress  exists,  and 
therefore,  neglect  the  transverse  normal  stress,  a^,  as  small  in  relation  to  the  other  stresses. 

2.  the  transverse  shear  strains,  exz,  Byz,  are  distributed  parabolically  through  the 
thickness  of  the  shell  and  the  transverse  shear  stresses,  a^,  Cyz,  satisfy  stress  free 
conditions  on  the  shell’s  top  and  bottom  surfaces. 

3.  the  shell  follows  linear  elastic  laminated  material  response  and  therefore  is 
restricted  to  small  strains. 

4.  the  shell  geometry  is  shallow;  Donnell  shallow  shell  strain-displacement 
relations  are  assumed. 

5.  the  shell  geometry  and  the  magnitude  of  the  loading  may  result  in  large 
deflections;  therefore,  assume  geometric  nonlinearity  in  the  von  Karman  sense. 

As  mentioned  in  the  introduction,  numerous  plate  and  shell  studies  have  shown 
that  due  to  the  relative  weakness  in  shear  of  current  high  stiffness-high  strength  laminated 
materials,  transverse  shear  deformation  can  be  a  significant  influence  even  for  relatively 
thin  geometries  as  compared  to  similar  isotropic  shells.  Consequently,  the  present 
approach  will  neglect  the  transverse  normal  stress,  a„  as  small  compared  to  the  shell’s 
inplane  stresses,  o^,  Oj,,  Gxy,  as  is  typical  in  plane  stress  assumptions.  However,  the 
approach  will  approximately  account  for  transverse  shear  deformation  by  retaining  the 
transverse  shear  stresses,  and  Oyz.  This  approach  is  consistent  with  shell  stress  order  of 
magnitude  studies  by  Koiter  (1967)  and  John  (1965)  for  isotropic  shells.  Additionally, 
exact  solutions  of  laminated  flat  plates  (Pagano,  1970)  and  cylindrical  shells  (Ren,  1987) 
in  cylindrical  bending  show  this  to  be  a  reasonable  approximation.  Furthermore,  since  Gz 
generally  is  of  order  h/R  times  the  bending  stress  and  Grz  und  are  of  order  hlr  or  his 
times  the  bending  stress,  Gz  is  negligible  compared  to  Grz  and  Gyz  when  r ,s  «R.  (Reddy, 
1985,.Palazotto  and  Dennis,  1992) 

Classical  Donnell  shallow  shell  assumptions  restrict  R/fc>500,  (Whimey,  1984). 
However,  Palazotto  and  Linneman  (1991)  show  that  for  shell  vibration  and  buckling,  the 
added  transverse  shear  degrees  of  freedom  allow  accuracy  to  approximately  R/h>50. 
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FIGURE  1.  Shell  geometry  of  radius  R,  thickness  h,  planform  rxs. 


KINEMATICS  AND  STRAIN 

Assuming  truncated  power  series  in  the  transverse  coordinate,  z,  for  the  shell 
displacements  and  satisfying  stress-free  conditions  on  the  shell’s  top  and  bottom  surfaces 
results  in  the  cubic  kinematics  of  Eqn  (1).  Eqn  (1)  gives  parabolic  approximations  for  the 
transverse  shear  strains  through  the  shell  thickness.  Kinematics  of  Eqn  (1)  and  similar 
forms  have  been  studied  by  several  investigators  (Levinson  1980,  Murthy  1981,  Reddy 
1984,  Bhimaraddi  1984,  Reddy  and  Liu  1985,  Soldatos  1987,  Dennis  and  Palazotto  1990) 
applying  them  to  both  laminated  flat  plate  and  shell  constractions  in  both  linear  and 
nonlinear,  buckling,  and  vibrational  analyses. 


Ui(x,y,z)  =  u+ z^i  +  ) 

(.x,y,  z)  =  V  -I-  z^2  +  ) 

U3ix,y)  =  w 


where,  Mi,  Uz,  M3  are  the  shell  displacements  in  the  x,  y,  z,  coordinate  directions 
respectively;  u,  v,  w  are  shell  midplane  displacements  (i.e.,  at  z=0);  and  ^1,  'F2  are 
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rotations  of  the  normals  to  the  shell  midplane  about  the  y,  x  axis  respectively,  and  commas 
denote  partial  differentiation  with  respect  to  the  coordinate  shown. 

Donnell  strain-displacement  relations  including  a  simple  geometric  nonlinearity  in 
the  von  Karman  sense  are  given  below  in  Eqn  (2)  (Brush  and  Almroth,  1975).  The 
nonlinearity  assumes  the  strains  and  rotations  squared  are  small  compared  to  one  and 
allow  deflections  of  magnitude  the  same  order  as  the  shell  thickness  (Chia,  1988).  The 
typical  stress  and  strain  shorthand  notation  defined  in  Eqn  (2)  is  adopted. 

W  1  2 

8,  =  £2  =  «2»y— 

£xy  =  £6  =  “l’>+«2..+>^»*^»y  (2; 

£yz  =£4  =  “2»j+*^3»)r 
£„  =  £5  =  “3».+«1». 

The  Donnell-von  Karman  strain-displacement  relations  of  Eqn  (2)  together  with  the  shell 
kinematics  of  Eqn  (1)  give  the  shell  strain-displacement  relations  in  Eqn  (3). 

£i=£;  +  zKi,  +  z'Ki3 
£2  =  £2  +  ^K21+z'K23 

£6  Kes  (3' 

£4  =  £>z'K42 

£5“£5‘*‘^  K52 

where  the  midplane  strains  are  given  by, 

£2  =  V,,-— +  iw,j, 

R 

£;  =  M,y+v„+w„w,, 

£4  =  '^2+^’y 

e;='Fi+w„ 


and  the  curvature  strains  are  given  by. 
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K21  =  ^2., 

Kei  “  '^i»>‘*'^2>» 


Ki3  =  ^/.(^Px+>»'.xc) 

K23  =  ^A('*'2.>+>^»yy) 

K  63  “  (^1  »>  '^'^2  >x  ■*‘2^>jit)r  ) 


K42  =  ^k,{%+W,y) 

Ks2  =  3A:*(yi+>v„) 


CONSTITUTIVE  RELATIONS 

Consider  shells  constructed  of  layers  of  unidirectional  filers  imbedded  in  a  matrix, 
i.e.,  transversely  isotropic  material.  (Jones,  1975)  The  three  dimensional  stress-strain 
relationships  for  a  single  lamina  of  transversely  isotropic  material  are  reduced  to  two 
dimensions  using  the  plane  stress  assumption.  By  neglecting  the  transverse  normal  stress 
(Oz=a3=0),  the  transverse  normal  strain  (£3 )  becomes  dependent  on  the  in-plane  strains 
(£1,82).  Eliminating  £3  then  gives  stress-strain  relationships  in  material  coordinates  where 
the  fibers  of  the  lamina  are  aligned  with  the  x  coordinate,  within  an  individual  ply  of  a 
laminated  shell  as  shown  in  Eqn  (4).  Shell  stractures  are  layered  with  the  material 
governed  by  Eqn  (4)  where  the  fibers  may  be  oriented  arbitrarily  within  the  plane  of  the 
shell,  i.e.,  the  x-y  plane.  Transforming  Eqn  (4)  from  material  coordinates  to  shell 
coordinates  then  gives  Eqn  (5)  for  the  nth  ply  of  the  laminate. 


Gi 

ra, 

Qn 

0 

0 

0  ‘ 

r  ^ 

81 

a2 

Q. 

Q22 

0 

0 

0 

82 

06 

^  = 

0 

0 

G66 

0 

0 

86 

04 

0 

0 

0 

0 

84 

05 
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0 

0 
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where  the  stress  notation  is  defined  analogously  to  the  strains,  see  Eqn  (2)  and  the 

are  elements  of  a  symmetric  array  of  the  reduced  stiffnesses  and  can  be  written  in  terms  of 
the  engineering  constants  as  given  below, 

Q44~  Gzi’  Q55~ 


2n  Qn  2i6  fe, 

On  Qu  K 


e«  a 
a 


where  the  ^  are  elements  of  a  symmetric  array  of  the  transformed  reduced  stiffnesses 

and  are  functions  of  the  reduced  stiffnesses  and  angle  of  orientation  of  the  fibers  of  the  nth 
ply.  (Jones, 1975) 

A  two  dimensional  shell  approach  will  “smear”  or  average  the  constitutive 
relations  of  each  ply  by  integrating  Eqn  (5)  through  the  shell  thickness.  First  substitute  the 
strains  of  Eqn  (3)  into  Eqn  (5),  see  Eqns  (6).  Using  averaged  load  quantities  defined  in 
Eqns  (7),  the  load-strain  relations  of  ^ns  (8)  result  where  in  Eqns  (8),  repeated  indices 
imply  summation,  =  1,2,6;  m,n  =  4,5;  and  the  elasticity  arrays  are  defined  in  Eqn  (9). 


a.  Qn  a 

Qn  a 

a 


£i  Kn  Ki3 

8^  +  2^K2i'  +  z''K23 
£6  .Kei.  Kes. 


o4"=  *2*1  a,  [fe; 

-oJ "  a,  lie; 


Ml  [M.l  [Pi]  %  [o. 

mJ  [mJ  [pJ  ioi. 


QnRs\-l\or^ 


Mi  =  5v£;  +  AKyi  +  FvKy3 
Pi  =  A£y‘*' AKyi  + i/i,Ky3 


Q„=A^tn  +  Dnu,Yi, 

Rm  =  Dmn2,\'^FmnYi, 


(8) 


VIRTUAL  WORK  AND  EQUILIBRIUM 


The  principle  of  virtual  displacements  for  the  shell  loaded  transversely  by  traction, 
q,  is  stated  as. 


+02  ^82+06  ^£6+04  Sel+Gs 


J  qhwdSl  —  0 


(10) 


where  Q„  is  the  shell’s  x-y  plane  within  the  nth  ply. 


Substituting  the  quantities  defined  in  Eqns  (3)  and  (7),  the  virtual  work  expression 
can  be  written  in  two  dimensions  as  is  shown  in  Eqn  (11). 


Jj  (M-  5e:+M.-  5k„+P.-  5k„+M-  hl+M,-  5k2.+/’=-  5k23+ 

Q 

+  Sk61  +  P6-  5k63  +  24  +  Sk42+  (11) 

Qs  Ses + Rs  Sk52 - =  0 

The  shell  equilibrium  equations  are  found  from  Eqn  (1 1)  by  substituting  the  strain- 
displacement  equations  of  (3)  and  then  integrating  by  parts.  Since  the  variation  of  each  of 
the  five  displacement  functions  is  arbitrary  and  independent,  their  coefficients  can  each 
individually  be  set  to  zero.  With  this  in  mind,  terms  were  grouped  for  each  virtual 
displacement  and  the  resulting  five  equations  are  given  in  Eqn  (12).  The  five  equilibrium 
equations  are  a  subset  of  Eqns  (12)  where  equilibrium  within  the  shell  domain  is  described 
by  each  double  integral  term  vanishing.  The  remaining  terms  in  each  equation  are  the 
boundary  conditions  for  the  shell  where  along  the  indicated  shell  edge,  either  displacement 
is  prescribed  or  the  force  term  is  zero.  The  given  form  of  Eqns  (12)  wiU  be  useful  for  the 
subsequent  Galerkin  solution.  The  underlined  terms  of  Eqn  (12c)  are  known  as 
“buoyancy”  terms  (Chia,  1980)  and  vanish  due  to  the  inplane  equilibrium  of  Eqns  (12a, b). 
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Excluding  the  nonlinear  terms  of  Eqn  (12),  i.e.,  those  where  a  force  quantity  multiplies  a 
displacement,  gives  equihbrium  equations  similar  to  those  of  other  investigators.  (Reddy 
and  Liu,  1985;  Linneman  and  Palazotto,  1991) 

II  (N,..  +  Ne.y'pi*  =  0  (12a) 

Q  «  0 

II  (V:,, + V6..)5v  dfl  - 1  A(,5v[rfj’  - 1  =  0  (12b) 


-\\(N.w.^*N,^w,+NjK*N^w,„*N2.,w,+2NeW,„+N,.w,+^ 

n 

]{N,w,+N,w,,-Pv,K-iPs.,k,+&  )&{<()+ 

0 

](Nyyo„*Nsw,.-P2.,t,-2P,..k,->-Q,  +^kX  )5H{*r  + 

0 

2t, P,5w[  i*]k,P, 5w,J^<fr+ 1 k,P^ Sw,!*  =  1198^.10 

ll{Mu*M.,*k,P,,,+k,P,,,-Q,-ik,R,'P'¥,da 

Q, 

-]  (Mi+ k,P^  'p{«’  -1  (M.+  k,p^'^  '?,[*  =  0 

0  0 

ii(My.,->-M,y.*kAy,*k,P,-Q,-3k,R,)S’v,da 

-|(M6+*.  P.)^'rldy-]{M2*k,  ft  )5T{<fc  =  0 


(12c) 


(12d) 


(12e) 
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III.  Galerkin  Solution 

The  Galerkin  technique  is  one  type  of  the  general  methods  of  weighted  residuals 
used  to  find  approximate  solutions  to  differential  equations  (Sokolrukoff,  1956).  In  this 
section,  the  Galerkin  technique  is  described  first  generally  and  then  the  technique  is 
applied  to  the  shell  equations  of  (12). 

Assume  a  differential  equation  of  the  form  shown  in  Eqn  (13) 

Lm“/=0  (13) 


where  L  is  &  differential  operator,  u  is  the  field  variable  (displacement  for  the  present 
study),  and  f  is  the  loading.  Next,  approximate  m  by  m  as  in  Eqn  (14) 

u=cJ^jyj  =  l,2X-N  (14) 

where  there  is  an  implied  summation  on  j,  the  c  .  are  unknown  constants,  and  the  ([)  . 

J  •' 

are  coordinate  functions  that  satisfy  the  boundary  conditions.  Substituting  the 
approximate  displacement  of  Eqn  (14)  into  (13)  then  gives, 

Lu-/=^  (15) 


In  Eqn  (15)  the  differential  equation  does  not  equal  zero;  but  instead  equals  some  error 
due  to  the  approximation  on  the  field  variable.  In  the  weighted  residual  techmque,  the 
constants  c  .  of  the  approximation  of  Eqn  (15),  are  found  such  that  the  error  in  Eqn  (15) 

is  forced  to  zero  in  an  average  sense  by  integrating  Eqn  (15)  multiplied  by  weighting 
functions  over  the  domain  as  in  Eqn  (16). 

J  (Lm-/)(1/  idD=0,  i=l,2,...N  (16) 

D 


The  case  where  the  weighting  functions,  \j/  j- ,  are  identical  to  the  coordinate  functions  of 

Eqn  (14)  is  the  Galerkin  method.  Eqn  (16)  represents  N  simultaneous  equations  in  the  N 
unknowns,  Cj.  Therefore,  a  single  differential  equation  is  approximated  to  desired 

accuracy  by  N  algebraic  equations. 

The  present  shallow  shell  approach  can  be  posed  in  the  operator  form  of  Eqn  (13) 
as  given  below  in  Eqns  (17)-(19).  In  the  equations,  the  force-displacement  relations  of 
Eqn  (8)  with  Eqn  (3)  are  substituted  into  the  equilibrium  equations  of  (12)  giving  five 
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coupled  partial  differential  equations  in  displacement  A  significant  simplification  resulte 
for  symmetrically  stacked  laminates  where  the  elasticity  arrays  of  Eqn  (9)  associated  with 
odd  powers  of  z  are  zero.  The  symbolic  math  package,  Mathcmaticu  (W^olfram  Research 
Inc.  1993),  performed  the  extensive  algebra. 

[L+H)d  =  f  (n: 

where  L  and  //  are  the  linear  and  nonlinear  operators,  respectively,  and  are  defined  as, 


jl2 


0  0 


Li2  L22  L23  0  0 

Ln  Loi  Lai  L34  Lis 

0  0  £ji4  Lu  Las 

0  0  L,,  L.S  Lss 


"0  0 

0  0 

//=  0  0 
0  0 


//o  0  O' 

//t,  0  0 

-H„  0  0 
0  0  0 


0  0  0  0  0 


and  the  elements  of  the  operators  are  defined  below, 

LlX  )“Al(  )>xy'b''^6(  )’yy 

Llli.  )“A6(  )>xi:'*'[A2  )»Jiy'^A6(  )»xy 


)“  A6(  )>X)r‘*‘A2(  )») 

Ln()  =  -^(  t 
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)=#(  )-[2A^-*-nD„K+KF„K^l  \„-{Ass*6D,,k,+9P„K%  I 
-[A^  +  6D^k,  +  9F^K^l  ),^+H22K''(  >>w 

■*^■^16^/1  (  ^xxxy'^^IS^h  (  ^xyyy'^^ee^h  (  )>x<3r)i 

L34(  )  =  -[Ass  +  6Dssk,+9Fssk,^l  ),.-[A,s  +  6D,sK+9F,,k,^l  >, 

+[F,^  +H,A"l  \:^+[F^K  >^+[3Fi,fc,  +3//A't  U 

+[FiA  +2F,,k,  +H,,K^+2H^k,^l  ),^ 

Lssi  )  =  -[AA5+^D,,k,+9F,,k,^l  ),.-[A^+6D^k,+9F^k,^l  >, 
+[F,,k,+H,,k,^l  ),^+[F22k,+H,,k,^l 

+[^12^A  ■*■  ^Fggfc^  + ](  )>*ir)r'^[3^26^*  ^  )»J5y 

L44(  )  =  4''5S  +9^'55^."](  )+[Ai  +2^. A  +«.!<:/](  >«. 

+[22^16 +2Hj5A:j^](  ^;,,+[D(6+2F66fc*  +  W(*^fc  ](  ),„ 

L«(  )=H;^«  +62>4st.  +9^'4i*.'](  )+[A.  +2F,A  k 

+[^12  ^66  ^I2^h  ^  ^66^h  ](  )>a:)' 

'^\P26^^^26^h'^  ^26^h  ](  ^yy 

Lssi  )  =  -[^44+6^44^*+9V/](  )+K+2FA+^66^*'](  >« 

■•■[2^26  +2i/26^/i^X  )>jy‘*’[^22  ■*■2^22^*  ’*'^22^*  ](  )>xy 

Hlsi  ')~^>x  Llli  Llli  ) 


Hlii  )  ^'•x  Llli  J^22(  ) 

//.,(  )  =  («..+|wJ)l.(  )+(v„-^+i«’']L.(  )+ 

{u,y+V,x+W,xW,y)L3i  LlsC  )-^^>yL23(  ) 


Lli  )~  Al(  )>»‘*‘^A6(  )»iy'*'A2(  )’»' 


T j-^i.  )~^12(  )»xx‘^^^26(  ^yy 

T j^(  )“^16(  )>xc'*'^^66(  )>jy'*’^26(  )>yy 


The  displacement  and  load  vector  are  given  by, 


(20) 


SIMPLY-SUPPORTED  BOUNDARY  CONDITIONS 

Consider  the  simply-supported  boundary  conditions  given  below  in  Eqn  (21). 

These  give  stress-free  membrane  conditions  normal  to  the  shell  plane  edges;  yet,  stress 
may  develop  tangentially  along  each  edge  due  to  the  restrained  tangential  displacement 

along  edges  x=0,  x=r: 

V  =  w  =  ^2  =  0 

along  edges  y=0,  y=s:  (21) 

u  =  w  =  'Pj  =  0 

In  the  Galerkin  approach,  the  displacements  of  the  shell  are  approximated  in  the 
form  of  Eqn  (14)  below  in  Eqn  (22),  The  coordinate  functions  must  satisfy  the  geometric 
boundary  conditions  of  Eqn  (21)  and  presently,  a  double  Fourier  series  is  assumed  for 
each  of  the  five. 


13 


where  the  coordinate  functions  are. 


,  nmx  .  nny 

=cos - sin - 

.  mux  my 

Xm„  =sin - cos - 

r  s 

.  rmvc  .  nicy 

=sin - sin - 

r  s 

nmx  .  my 

=cos - sin - 

r  s 

.  nmx  my 
Tl™,  =siit - cos - 

» mn  _ 


Substituting  Eqn  (22)  into  Eqn  (17)  plus  reintroducing  die  nonzero  boundary 
condition  terms  of  Eqn  (12),  gives  the  Galerkin  equations  of  the  form  of  Eqn  (16)  in  Eqns 
(23).  Since  admissible  functions  are  assumed  for  the  displacements,  i.e.,  those  that  only 
satisfy  the  geometric  boundary  conditions,  the  single  integral  terms  are  the  remaining 
nonzero  boundary  conditions  and  must  be  retained  as  is  shown. 


jl[Lii(«)  +  Li2(v)  +  Li3H  +  /fi3M>w^y  + 

0  0 

I  Ni  (0,  >’>l>  „  (0,  y)dy-]  N I  (r,  y)(l)  (r,  y)dy  =  0 


(23a) 


1 1  [Li2  (w)  +  L22  (v)  +  L23  i^)  +  if  23  (>*')]  X  p,dxdy  + 

0  0 

r  r 

J  jVz  U.0)x  „  (X,0)  d!r  - 1  ‘^)X  pq  (x,  s)dx  =  0 


(23b) 


jl[Ll3(M)  +  L23(v)  +  L33H  +  L34(^l)+L35(^2)-ii33W]a,,^^3'  + 

0  0 

0  0 

r  r  r  s 

K\ {x,s)dx-k^\ p^ixfiyx  ixfi)dx  =  j  J  qa^,^  dxdy 


(23c) 
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0  0 

jMi(0,}')Yp,(0,}')rf3'-jMi(^>’)Yp,(^3’)^>’+  (23d) 

0  0 

^*1  Pi(0,yyr  p,iO,y)dy-k,j  Pi(r,yyr  p^ir,y)dy  =  0 
0  0 


0  0 

I  M2(^’0)n,,(x,o)dx  -  j 

0  0 

k]  P2ix,0)np,(x,0)dx-kjj  p,(x,s)(\^^ix,s)dx  =  0 
0  0 

Eqns  (22)  are  next  substituted  into  Eqns  (23)  giving  a  system  of  5MN 
simultaneous  nonlinear  algebraic  equations  in  the  unknown  constants  (imn>  bmn>  Cmm  dmn> 
Cr^.  In  getting  to  the  simultaneous  equations,  many  trigonometric  integrals  must  be 
evaluated.  The  general  forms  are  shown  in  Appendix  A.  The  final  form  of  the  Galerkin 
equations  given  as  functions  of  displacement  is  found  in  Appendix  B.  Appendix  C  is  a 
listing  of  the  Fortran  coding.  In  addition,  sample  input  and  output  listings  are  included. 


LOADING 

Various  transverse  loadings  can  be  considered  by  first  expanding  the  traction  q  of 
the  right  hand  side  of  Eqn  (23c)  into  a  Fourier  series  as  shown  in  Eqn  (24). 


g  =  gix,y)  =  an-— sm— 

—1  f  s 


Following  Timoshenko  and  Woinowsky-Kreiger  (1959),  any  coefficient  of  the  series  can 
be  found  from: 


4  f  f  .  .  .  iTCC  .  jKy,  , 

Pj  =  —  J  J  ^(j:,y)sin — sm - dxdy 


The  following  transverse  loadings  are  then  derived  from  Eqns  (24)  and  (25), 
Uniform  loading  of  magnitude  q,: 
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(26) 


(27) 


(28) 
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rV.  Linear  Plate  and  Shell  Solutions 

As  a  first  step  in  verifying  the  solution  procedure,  several  linear  cases  are 
presented  and  compared  to  published  results.  Linear  solutions  assume  proportional  load 
and  deflection  and  therefore  the  terms  of  the  operator  H  in  Eqn  (17)  are  not  included.  For 
many  of  the  cases  examined,  biaxial  symmetry  is  present  and  consequently  only  the  odd 
terms  of  the  assumed  displacement  functions  of  Eqn  (22),  i.e.,  m,n  =  l,3,5...etc., 
contribute  to  the  solutions,  (^uasi-isotropic  plates  and  shells,  such  as  the  [0/±45/90]s  and 
[-60/0/60]s  laminates  examined  in  this  section,  have  nonzero  Die  and  1)26  terms  for 
example  that  destroy  the  symmetry  and  in  that  case,  both  even  and  odd  displacement 
functions  are  required.  Unique  results  are  given  for  a  [-60/0/60]s  shell  laminate  using  tire 
recently  proposed  Batdorf-Stein  parameter. 

ISOTROPIC  PLATE 

Consider  a  thin  isotropic  plate  subjected  to  two  different  transverse  loadings, 
uniform  pressure  q,  and  a  single  point  load  Po  applied  at  the  center  of  the  plate,  see  Eqns 
(26)  and  (28).  For  a  thin  geometry,  i.e.  for  this  case,  the  thickness  ratio  S=.s//i=100, 
transverse  shear  deformation  is  negligible  and  comparison  to  the  classical  thin  plate 
response  should  be  close.  Table  1  gives  the  center  plate  transverse  displacement  solution 
for  an  increasing  number  of  terms  in  the  assumed  solution  of  Eqn  (22).  As 
mentioned,  for  uniform  and  center  plate  point  loads  and  isotropic  inaterial  behavior,  only 
the  symmetric  coefficient  functions,  i.e.,  odd  m,n  in  Eqn  (22)  contribute  to  the  solution. 
Consequently,  m  is  1,3,5,...M  and  n  is  1,3,5,..  JV.  Both  of  the  uniform  load  cases  converge 
to  within  0.1%  of  the  classical  result  taken  from  Pilkey  and  Chang  (1978)  after  only  a  few 
terms  in  the  series.  The  point  load  case  requires  several  additional  terms. 

CROSS-PLY  LAMINATED  PLATE 

Results  are  shown  for  square  cross-ply  laminates,  [0/90/0]  and  [0/90]®.  The 
[0/90/0]  laminate  is  subjected  to  a  uniform  load  of  Eqn  (26)  and  the  nondimensionalized 
center  plate  deflection  of  Eqn  (29a)  requires  several  terms  for  convergence,  see  Table  2. 
The  increasing  flexibility  due  to  transverse  shear  deformation  is  evident  in  Table  2  for  the 
smaller  thickness  ratios,  S,  For  and  odd  the  present  Galerkin  solution 

exactly  matches  the  results  due  to  Reddy  (1984). 

The  [0/90]s  laminate  is  subjected  to  the  sinusoidal  pressure  loading  of  Eqn  (27) 
and  one  term  due  to  this  special  loading  case,  gives  the  exact  solution.  Table  3 

shows  the  nondimensionalized  deflection  and  stress  of  Eqns  (29).  All  results  of  Table  3 
are  identical  to  those  given  by  Reddy  (1984).  Transverse  shear  deformation  is  again 
evident  in  the  thick  plates. 
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TABLE  1.  Center  plate  deflection  for  an  isotropic  thin  plate.  Deflection  is 
nondimensionalized  by  lOQDtqa^  for  uniform  loads  and  by  lOOD/Pfl^  for  the  point  load 
case,  where  D=E/i^/i2(i '0=0.3;  Square:  r=j=fl;  rectangular:  r=2s=2a. 


square 

rectangular 

square 

M,N 

uniform 

uniform 

point 

load 

load 

load 

1 

0.41607 

1.06514 

1.02720 

5 

0.40637 

1.01387 

1.14332 

9 

0.40624 

1.01297 

1.15485 

13 

0.40624 

1.01288 

1.15815 

29 

1.01287 

1.15897 

classical  (MJ^  =99 

0.40624 

1.01287 

1.16002 

in  Pilkev,  1978) 

* 

where  in  Eqns  (29),  from  Figure  1,  r  and  j  represent  the  shell  planform,  /i  is  the  shell 
thickness,  and  q  is  defined  in  Eqns  (24)-(28);  the  quantities  are  evaluated  at  the  plate 
coordinates  given  in  parentheses. 


18 


TABLE  2.  Nondimensionalized  center  plate  deflection  for  [0/90/0]  laminate 
subjected  to  uniform  load;  EilE2=25,  Gi2/£2=Gi3/£^2=0.5,  G23/E2=0.2,  'Ui2=0.25. 


S=s/h 

MM 

1  5  13  21 

Reddy 

(1984) 

100 

0.7039  0.6709  0.6705  0.6705 

0.6705 

50 

0.7182  0.6843  0.6838  0.6838 

0.6838 

20 

0.8173  0.7769  0.7760  0.7760 

0.7760 

10 

1.1551  1.0921  1.0901  1.0900 

1.0900 

4 

3.1156  2.9167  2.9094  2.9091 

2.9091 

2 

8.3141  7.7823  7.7665  7.7661 

7.7661 

TABLE  3.  Nondimensionalized  center  plate  deflection  and  stress  for  [0/90]s  laminate 
subjected  to  sinusoidal  pressure;  EilE2—25,  Gi2/E2=Gi2lE2=0.5,  G2ilE2=0.2,  \)i2=0.25; 
MJ^=L 


S=s/h 

w 

Cl 

02 

04 

05 

100 

0.4343 

0.5387 

0.2708 

0.1117 

0.2897 

0.5060 

0.5393 

0.3043 

0.1234 

0.2825 

0.7147 

0.5456 

0.3888 

0.1531 

0.2640 

4 

1.8937 

0.6651 

0.6322 

0.2389 

0.2064 

CROSS-PLY  CYLINDRICAL  SHELL 


Similar  to  the  flat  plate  case,  a  [0/90]s  cross-ply  cylindrical  shell  subjected  to 
transverse  sinusoidal  pressure  is  solved  exactly  by  one  term  in  the  assumed 

displacement  of  Eqns  (22).  The  present  solution  is  compared  to  a  solution  given  by  Reddy 
(1982)  who  assumed  the  lower-order  Reissner-Mindlin  transverse  shear  deformation 
theory.  Table  4  compares  the  two  approaches  for  a  thick  shell  geometry,  S=s/h=l0  and 
two  curvature  ratios,  R/h=10  or  100.  The  present  Galerkin  approach  gives  a  more  flexible 
response  for  both  R/h  ratios  as  expected  due  to  its  more  accurate  transverse  shear 
representation.  Also  shown  is  the  solution  from  a  mesh  of  shear  deformable  finite 
elements  by  Palazotto  and  Dennis  (1992).  The  finite  element  solution  is  also  based  on  the 
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parabolic  transverse  shear  distribution  through  the  shell  thickness.  However,  although  the 
Galerldn  and  Reddy  results  are  based  on  Donnell  shallow  shell  theory,  the  finite  element 
solution  can  model  deep  shell  behavior.  Even  for  the  deeper  shell  (R/h=10),  the  present 
Donnell  Galerldn  equations  give  fairly  accurate  results  as  compared  to  the  deep  shell  finite 
element  solution. 

The  shell  Batdorf-Stein  parameter  Z,  defined  in  Eqn  (30)  for  balanced  S3mimetnc 
laminates,  is  useful  in  shell  studies  where  the  response  of  different  geometries  are 
compared  (Nemeth,  1994). 


i4.li  i422  i4: 


1/2 


(30) 


where  s  is  the  circumferential  length  of  the  shell.  Ay  and  Dy  are  defined  Eqn  (9). 

For  comparison  of  like  material  unidirectional  [0]  and  [90]  shell  laminates,  Z  reduces  to 
Eqn  (31).  Eqn  (31)  has  been  used  in  nondimensional  buckling  and  nonlinear  laminated 
shell  studies  by  Dennis  et  al  (1993,1994).  ShaUow  shells  are  indicated  by  small  vdues  of  Z 
and  the  flat  plate  limiting  case  has  Z=0.  The  values  of  Z  for  the  shells  of  this  section  are 
shown  in  the  Tables. 


Z=^[l (31 

Convergence  characteristics  of  the  present  Galerldn  solution  for  a  shell  geometry 
are  seen  by  subjecting  the  [0/90]*  shell  laminate  to  a  uniform  transverse  pressure.  Results 
are  shown  in  Table  5  for  odd  m,n  terms  only  and  are  compared  to  an  unpublished  finite 
element  solution. 

QUASI-ISOTROPIC  PLATES  AND  SHELLS 


Two  quasi-isotropic  geometries  are  subjected  to  a  uniform  load.  The  Galerldn 
solution  is  compared  to  the  published  results  of  Phan  and  Reddy  (1985)  in  Table  6  for  a 
[0At45/90]s  flat  plate  laminate  using  the  nondimensionalization  of  Eqn  (29).  The  Galerldn 
solution  is  based  on  including  terms  of  both  the  even  and  odd  /w,/i  of  Eqn  (22)  as  the 
nonzero  Die  and  D26  stiffnesses  due  to  the  ±45"  plies  destroy  the  biaxial  symmetry  that 
exists  in  the  previous  cases.  However,  the  Phan  results  assumed  symmetry  in  their  study 
and  this  accounts  for  the  differences  with  the  present  where  for  every  thickness  ratio,  the 
Galerldn  solution  gives  a  stiffer  response.  An  independent  finite  element  analysis  based  on 
meshes  of  the  28  degree  of  freedom  flat  plate  element  developed  by  Palazotto  and  Deniiis 
(1992)  confirmed  this  where  center  plate  deflections  of  Table  6  were  also  found  both  with 
and  without  the  symmetry  assumed. 

Secondly,  deflection  and  stress  are  shown  in  Table  7  for  a  [-60/0/60]s  quasi- 
isotropic  shell  laminate  for  several  geometries.  Again,  the  even  m,n  terms  in  the  assumed 
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displacement  are  required  due  to  the  nonzero  D\6  and  D26  terms.  For  a  given  value  of  the 
Batdorf-Stein  parameter  Z,  the  thicker  shells  represented  by  smaller  values  of  S  undergo  a 
more  flexible  response.  Transverse  shear  flexibility  accounts  for  the  larger  deflections. 

For  a  given  thickness  ratio  sih,  as  the  curvature  is  increased  as  indicated  by  increasing  Z, 
both  the  center  shell  deflection  and  circumferential  stress  decrease  due  to  more  membrane 
action  and  less  bending.  These  results  are  also  consistent  with  the  response  given  in  Table 
4.  For  many  of  the  geometries  with  larger  Z,  the  thicker  shell  responses  (smaller  S)  are 
not  shown  as  Rlh  becomes  too  small. 

TABLE  4.  Nondimensionalized  center  shell  deflection  and  stress  for  cross-ply  cylindrical 
shell  subjected  to  sinusoidal  transverse  pressure;  EilEz=Q.5,  GijJE2-GiilEi=Q.5, 
G23lE2=0.2,  -012=0.25;  S=s/h=r/h=10;  MJ^=h 


w 

Z=  1.177 
R/h=m 

Oi 

02 

w 

Z=  11.77 
R/h=10 

Oi 

O2 

Galerkin 

0.11231 

0.5501 

0.3938 

0.53645 

0.4570 

0.3393 

Reddy (1982) 

0.6609 

0.4998 

0.3637 

0.5006 

0.4233 

0.3190 

Finite 

Element* 

0.71354 

0.5520 

0.3923 

0.55665 

0.4751 

0.3075 

elements  developed  by  Dennis  and  Palazotto  (1992). 


TABLE  5.  Nondimensionalized  center  shell  deflection  and  stress  for  cross-ply  cylindrical 
shell  subjected  to  uniform  transverse  pressure;  EilE2=25,  Gi2/Ei=Gi3lE2=0.5,  G23/£^2=0.2, 
012=0.25;  Z=11.77,/?/;i=1000,  s/h=100. 


1 

M,N 

5 

9 

13 

Finite 

Element* 

w 

0.5858 

0.5663 

0.5661 

0.5661 

0.5662 

Oi 

0.7785 

0,7302 

0.7286 

0.7283 

0.7293 

02 

0.4172 

0.3378  - 

0.3339 

0.3332 

0.3324 

♦Finite  element  solution  based  an  8x8  quarter  model  mesh  of  36  degree  of  freedom  shell 
elements  developed  by  Dennis  and  Palazotto  (1992). 
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TABLE  6.  Nondimensionalized  center  plate  deflection  for  [0/±45/90]s  laminate  subjected 
to  uniform  load;  £i/E2=40,  Gy^JEi^GnlEi^QA  023/1^2=0.5, 1)12=0.25;  m,n  =  1,2, 3, ...MA/ 


S=s/h 

m:n 

1  3  9  15 

♦Phan 

(1985) 

100 

0.3380  0.3354  0.3372  0.3377 

0.3769 

50 

0.3444  0.3412  0.3431  0.3437 

0.3859 

20 

0.3893  0.3818  0.3847  0.3852 

0.4336 

10 

0.5486  0.5259  0.5316  0.5321 

0.5904 

4 

1.6275  1.5070  1.5283  1.5289 

1.6340 

♦Assumed  plate  biaxial  symmetry. 


TABLE  7.  Nondimensionalized  center  plate  and  shell  deflection  and  stress  for  [-60/0/60]  s 
laminate  subjected  to  uniform  load,  Ei/£2=40,  Gi2JEz=GnlE2=0.6,  023/^2=0.5,  Di2=0.25, 
m,n  =  1,2,3,...21.  Stress  is  compressive  except  for  asteriked  value. 


Z 

S=slh  R/h  w  0  2  (z—hJ2) 

0 

100  -  0.3557  0.5149 

20  -  0.4006  0.5341 

10  -  0.5307  0.5072 

-5  -  1.0142  0.5949 

1.9384 

100  5000  0.3461  0.4657 

20  200  0.3878  0.4908 

10  50  0.5042  0.5337 

12.5  0.9191  0.6223 

4.846 

100  2000  0.2948  0.4426 

20  80  0.3238  0.4583 

10  20  0.3987  0.4761 

5  5  -  .  - 

9.692 

100  1000  0.1916  0.3298 

20  40  0.2019  0.3272 

10  10  0.2245  0.3063 

5  2.5  -  -  . 

24.230 

48.461 

96.922 

100  400  0.05199  0.1062 

100  200  0.01192  0.01800 

inn  100  0.0008239  ^0.01061 
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ANGLE-PLY  PLATE 

A  thick  plate  consisting  of  a  single  45°  ply  is  subjected  to  both  a  uniform  and 
sinusoidal  transverse  pressure  and  nondimensional  center  plate  deflection  and  stress  are 
shown  in  Table  8.  This  laminate  has  nonzero  Aie  and  A26  membrane  stiffnesses  and  its 
effect  apparently  slows  convergence  considerably  compared  to  other  results  of  this  section 
where  the  solution  including  24  terms,  both  even  and  odd,  is  yet  unconverged. 


TABLE  8.  Nondimensionalized  center  plate  deflection  and  bending  stress  {z=±hl2)  for 
angle-ply  [45]  Gi2/£2=Gi3/£2=0.5,-  G23/E2=0.2,  \)i2=0.25;  S=a//i=10,  a=r=s. 


load 

3 

6 

M,A 

12 

24 

uniform 

w 

0.6632 

0.7073 

0.7569 

0.7877 

Cl 

0.4609 

0.3621 

0.3885 

0.4175 

sinusoidal 

w 

0.4257 

0.4499 

0.4794 

0.4969 

Oi 

0.3198 

0.2552 

0.2749 

0.2909 
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V.  Nonlinear  Plate  and  Shell  Solutions 

Geometric  nonlinearity  becomes  important  for  plate  and  shell  cases  where  the 
transverse  displacement  is  no  longer  small  compared  to  the  thickness.  The  assumed  von 
Karman  nonlinearity  of  the  in-plane  strains  of  Eqn  (3)  is  valid  when  the  deflection  is  of  the 
same  order  of  magnitude  as  the  thickness  of  the  shell.  Consequently,  the  nonlinear 
solutions  presented  here  will  give  the  plate  and  sheU  response  up  to  several  shell 
thicknesses  of  deflection. 

An  iterative  algorithm  gives  the  solution  to  the  nonlinear  equilibrium  equations  of 
Appendix  B,  see  Eqns  (B1)-(B8).  The  algorithm  uses  direct  iteration  with  residual  forces 
determining  convergence.  Normally,  if  a  norm  of  the  displacement  vector  changed  by  less 
than  0.1%  from  the  previous  iteration,  convergence  was  assumed. 

Including  the  nonlinear  terms  of  Eqns  (B6)-(B8)  requires  numerous  calculations 
for  solution  of  the  shell  equilibrium  path.  To  gain  an  appreciation  for  the  additional 
computation  burden  of  the  nonlinear  solution,  compare  the  number  of  passes  through 
loops  of  computer  coding  for  the  linear  versus  the  nonlinear  solution.  Although  the 
number  of  matiiematical  operations  is  different  for  a  loop  within  the  linear  versus  the 
nonlinear  coding,  counting  the  number  of  passes  through  loops  will  stiU  illustrate  the 
significant  additional  math  required.  Let  n  equal  the  number  of  terms  in  one  of  the 
summations  in  the  assumed  displacement  of  Eqn  (22).  Calculating  aU  stiffnesses  of  the 
linear  coefficient  matrix  of  Eqns  (B1)-(B5)  requires  n'^  passes  through  algorithm  loops. 
The  nonlinear  coefficient  matrix  is  a  sum  of  the  linear  stiffnesses  that  are  calculated  only 
one  time  and  the  stiffnesses  of  Eqns  (B6)-(B8).  Calculation  of  the  stiffnesses  of  Eqns 
(B6)-(B8)  requires  n®  loops/or  each  iteration  of  the  nonlinear  algorithm.  Several 
iterations  are  typically  required  for  convergence  for  each  increment  of  load  that  represents 
a  single  point  on  the  nonlinear  equilibrium  path.  Since  accurate  solutions  were  typically 
found  using  nine  terms  in  both  x  and  y  coordinate  directions  (for  81  terms)  for  each  of  the 
assumed  displacements  of  Eqn  (22),  the  nonlinear  solution  algorithm  consumes 
significantly  more  computing  time  over  the  linear  solution.  In  cases  where  biaxial 
symmetry  exists,  the  solutions  given  resulted  from  m,n  =1,3,5,.. .17  in  the  series  of  Eqn 
(22). 

The  nonlinear  Galerldn  solution  technique  is  verified  by  calculating  the  response  of 
an  isotropic  plate  subjected  to  uniform  pressure.  The  classical  solution  was  published  by 
S.  Levy  in  1942  and  is  a  typical  nonlinear  plate  comparison.  Next,  plate  and  cylindrical 
shell  solutions  are  presented  for  unidirectional  and  cross-ply  laminates.  The  effect  of 
transverse  shear  deformation  in  the  presence  of  geometric  nonlinearity  is  discussed. 


ISOTROPIC  PLATE 

A  square  isotropic  plate  is  subjected  to  uniform  transverse  pressure  and  the 
response  is  plotted  in  Figure  2.  The  present  response  of  Figure  2  resulted  from 
m,n=l,3,5,...17  in  the  assumed  displacement  of  Eqn  (22).  For  increasing  load,  membrane 
resistance  results  in  a  stiffer  plate  response  in  the  nonlinear  solution  compared  to  the  linear 
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response  where  no  membrane-bending  coupling  is  modelled.  Presently,  the  nonlinear 
response  diverges  from  the  linear  solution  for  deflections  larger  than  approximately  one 
third  of  the  plate  thickness.  Although  the  response  shown  compares  well  to  the  published 
results  of  Levy  (1942)  for  the  load  range  given,  at  higher  loads  the  solution  begins  to 
diverge  from  the  published  and  increasingly  more  terms  in  the  assumed  displacement  are 
required  for  accuracy. 
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FIGURE  2.  Center  plate  deflection  versus  load  for  a  square  {r=s=a)  isotropic  plate 
subjected  to  uniform  pressure;  '0=0.316. 


UNIDIRECTIONAL  PLATE  LAMINATE 

Next  consider  a  plate  constructed  with  fibers  aligned  the  x  coordinate,  i.e.,  a  [0] 
laminate,  subjected  to  sinusoidal  transverse  loading  of  in  Eqn  (27).  The  response  is  shown 
in  Figure  3  for  two  plate  thickness  ratios  {S=alh,  r=j=a).  The  linear  response  for  both 
plates  is  also  shown  in  the  figure.  Due  to  the  simple  sinusoidal  loading,  a  nearly 
converged  response  results  by  including  only  a  few  terms  in  the  assumed  displacement.  In 
the  previous  uniform  load  case,  both  the  load  and  the  response  must  be  approximated  by 
the  trigonometric  series.  In  this  case,  the  load  is  modelled  exactly  by  only  one  term  and 
hence  easier  convergence  in  the  displacement  response  is  seen. 
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The  thick  plate  undergoes  a  more  flexible  response  compared  to  the  thin  plate  due 
to  increased  transverse  shear  flexibility.  However,  transverse  shear  deformation  in  the 
nonlinear  analysis  is  less  than  in  the  linear  analysis  since  in  the  nonlinear,  membrane 
stiffemng  is  a  competing  influence.  For  this  case,  whereas  the  linear  response  is  15.3% 
more  flexible  in  the  thick  plate  compared  to  the  thin  plate  for  all  load  levels,  the  nonlinear 
is  only  13.4%  greater  for  nondimensionalized  load  of  250, 11.1%  greater  for  a  load  of 
500,  and  only  10.6%  greater  for  a  nondimensionalized  load  of  750. 

Figure  4  shows  the  center  plate  normal  stress,  CTi,  for  the  thin  [0]  laminate,  linear 
versus  nonlinear.  The  stress  calculated  at  the  plate  center  for  the  thick  plate  is 
approximately  the  same  as  that  shown  in  Figure  4,  the  thin  response  being  slightly  higher 
than  the  thick  plate  response. 
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FIGURE  3.  Center  plate  deflection  versus  load  for  a  square  (r=s=a)  unidirectional 
laminate  subjected  to  sinusoidal  transverse  pressure;  EiIE^fAQ,  Gi2/E2=Gi3/E2=0.6, 
G23/E2=0.5,  1)12=0.25;  M^=Y1  (odd  m,n  only). 


CROSS-PLY  PLATE  LAMINATE 

The  response  of  a  [0/90]  §  laminate  under  uniform  transverse  pressure  is  given  for 
two  plate  thickness  ratios  in  Figure  5.  As  was  true  for  the  unidirectional  plate,  the  thicker 
plates  undergo  a  more  flexible  response  compared  to  the  thin  plate  due  to  transverse  shear 
deformation.  The  flexibility  due  to  transverse  shear  is  less  in  the  nonlinear  response 
compared  to  the  linear  response  as  membrane  stiffening  is  a  competing  influence  similar  to 
that  seen  in  the  unidirectional  case. 
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FIGURE  4.  Center  plate  normal  stress  in  unidirectional  square  (r=s=a)  laminate 
subjected  to  sinusoidal  transverse  pressure;  Ei/E2=40,  Gi2/E2=Gi3/E2=0.6,  G23/jE^2=0.5, 
1)12=0.25;  5=a//i=100;  M//=17  (odd  m,n  only).  Stress  extreme  tensile,  z=hl2. 
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FIGURE  5.  Center  plate  deflection  versus  load  for  a  square  (r=s-a)  cross-ply  plate 
subjected  to  uniform  pressure;  £i/£2=40,  Gi2/E2=Gn/E2=0.6,  G23/F2=0.5,  Di2=0.25; 
MJ^=n  (odd  m,n  only). 


27 


UNIDIRECTIONAL  SHELL  LAMINATE 

Unidirectional  [0]  and  [90]  laminates  are  subjected  to  sinusoidal  transverse 
pressure.  Figure  6  gives  the  center  shell  response  of  a  [90]  laminate  for  the  Batdorf-Stein 
parameter  of  Eqn  (30),  Z=9.992.  Initially,  the  shell  increasingly  deflects  with  load  due  to 
the  flexibility  gained  through  circumferential  membrane  compression.  However,  the 
enlarging  longitudinal  tension  field  eventually  overcomes  the  flexibility  gained  through 
membrane  compression  and  the  shell  begins  to  stiffen  with  load.  For  all  loads,  the  thicker 
shell  (,R/h=40,  alh=2G)  is  more  flexible.  The  response  of  a  [0]  laminate  was  also 
calculated  and  for  this  value  of  Z,  the  response  was  nearly  identical  to  the  [90]. 
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FIGURE  6.  Center  shell  deflection  versus  load  for  a  square  ir=s=a)  unidirectional  shell 
laminate  subjected  to  sinusiodal  transverse  pressure;  £i/£2=40,  Gi2/E2=Gi3/E2=0.6, 
G23/E2=0.5,  1)12=0.25;  (odd  m,n  only). 

CROSS-PLY  CYLINDRICAL  SHELL  LAMINATE 

A  thin  cross-ply  shell  is  subjected  to  sinusoidal  transverse  pressure  and  compared 
to  a  finite  element  solution  in  Figure  7.  Due  to  membrane  compression,  the  shell  initially 
exhibits  more  flexible  behavior  as  the  load  is  increased,  i.e.,  the  shell  softens  with  load. 
Eventually,  significant  deflection  results  when  the  load  is  increased  only  a  small  amount  as 
given  by  the  shallow  slope  in  the  load-deflection  curve.  The  longitudinal  bending  tension 
field  on  the  bottom  of  the  shell  enlarges  until  its  stiffening  influence  overtakes  the 
flexibility  caused  by  membrane  compression  and  the  shell  stiffens  for  the  highest  loads 
shown  in  the  figure.  The  Galerkin  solution  follows  the  finite  element  result  well  except  the 
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latter  reaches  a  slightly  higher  load  before  the  very  shallow  slope  in  the  load-deflection 
curve  begins. 

Figure  8  shows  the  nonlinear  response  of  two  cross-ply  shells  that  have  a  Batdorf- 
Stein  parameter  equal  to  1 1.945.  One  shell  geometry  is  thin  and  shallow  and  should 
experience  insignificant  transverse  shear  deformation.  The  second  shell  is  a  deeper, 
thicker  shell  where  the  more  flexible  response  shown  in  the  figure  is  probably  due  to 
transverse  shear  deformation.  Although  the  geometry  of  the  shell  (R,  h,  t,  5)  is  the  same 
as  for  the  unidirectional  shell  of  Figure  6,  the  Batdorf-Stein  parameter  of  Eqn  (30)  is 
larger  (9.992  versus  1 1.945)  due  to  the  ply  stacking.  Interestingly  however,  boA  thin  and 
thick  shell  responses  are  nearly  identical  for  the  unidirectional  and  cross-ply  laminates;  and 
furthermore,  have  the  same  value  Z  if  Eqn  (31)  were  used  where  ply  stacking  is  not 
considered  in  the  parameter  definition. 
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FIGURE  7.  Center  shell  deflection  versus  load  for  a  square  (j-=s=a)  cross-ply  shell 
laminate  subjected  to  sinusiodal  transverse  pressure;  £i/E2=40,  GvJE2^iilE%=Q£, 
G23/£2=0.5,  1)12=0.25;  M^=\l  (odd  m,n  only). 
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FIGURE  8.  Center  shell  deflection  versus  load  for  square  (r=s=d)  cross-ply  shell 
laminates  subjected  to  sinusiodal  transverse  pressure;  Ei/£2=40,  GuJEz^GijIEz^^.S, 
G23lE2=0.5, 0)12=0.25;  MJ^=n  (odd  m,n  only). 
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VI.  Conclusions 


A  laminated  shallow  shell  approach  that  includes  von  Karman  geometric 
nonlinearity  and  parabolic  transverse  shear  deformation  has  been  posed  in  differential 
operator  form.  Trigonometric  functions  are  then  assumed  for  each  of  the  five  shell 
displacement  degrees  of  fi’eedom  for  the  subsequent  nonlinear  galerkin  solution.  The 
galerkin  nonlinear  solution  is  computationally  intensive.  The  response  of  several  laminate 
geometries  subjected  to  transverse  loadings  are  found.  Thicker  plates  and  shells  generally 
exhibit  more  flexible  response  compared  to  thinner  geometries  in  both  linear  and  nonlinear 
analyses.  The  nondimensional  shell  response  is  examined  by  using  the  Batdorf-Stein  shell 
parameter  for  laminated  constructions.  Quasi-isotropic  shallow  shells  undergo  significant 
transverse  shear  flexibility  in  the  thicker  geometries  as  given  by  the  nondimensional  shell 
crown  deflection.  However,  the  nondimensional  crown  deflection  in  the  deeper  shell 
response  is  not  much  influenced  by  shell  tWckness.  For  flat  plates,  geometric  nonlinearity 
lessens  the  influence  of  transverse  shear  flexibility  when  compared  to  linear  solutions  due 
to  membrane  stretching  resistance. 
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APPENDIX  A 

TRIGONOMETRIC  INTEGRALS 

Evaluation  of  the  trigonometric  integrals  resulting  from  the  assumed  displacement  of  Eqn 
(22)  substituted  into  Eqn  (23)  makes  use  of  the  trigonometric  identities  of  Eqn(Al). 

cos  A  cos  [cos(A  -B)  +  cos(A  +  B)] 

sin  A  sin  B  = — [cos(A  —  B)  —  cos(A  +  B)] 

2 

sin  A  cos  B  = — [sin(A  -  B)  +  sin(A  +  B)] 

2  .  (Al) 

Linear  stiffness  terms  require  evaluation  of  the  three  integrals  of  Eqns  (A2)-(A4): 

Note:  i,j,  k,  /,  are  integers  and 


r 

J  cosix*  cos  jx'  dx  = 


2  if  i  =  J 


0  if  I  9^  j 


r 

Jsinix'sin  jx'  dx  = 


2  if  i  =  J 


0  if  i^j 


r 

J  sin  Cc"  cos  dx  = 


0  if  i  +  J  =  even 


Nonlinear  stiffiiess  terms  require  evaluation  of  the  four  integrals  of  Eqns  (A5)-(A8). 

r . ,  .  ,  .  .  r  1  1  1 1 

sin  ix  sin  /X  sin  kx  ox  =  —“(—+ - ^ 
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if  (i-j)-ik-l)  =  0 
if  otherwise 
if  (i-j)+(k-l)  =  0 
if  otherwise 
if  ii-j)-(k+l)  =  0 
if  otherwise 
if  (i-j)+ik  +  l)  =  0 
if  otherwise 
if  (i+7)-(A:-/)  =  0 
if  otherwise 
if  (/+7)+(A:-/)  =  0 

if  otherwise 

if  (i  +  7)-(*+/)  =  0 
if  otherwise 


APPENDIX  B 
GALERKIN  EQUATIONS 


Eqns  (23)  are  given  below  in  terms  of  constants  Oij,  by,  Cij,  dij,  and  Cij.  Five 
equilibrium  equations,  (B1)-(B5),  result  for  every  p,q  of  l^n  (23).  The  summations 
extend  to  M,N  see  Eqn  (22). 

Linear  terms 

Note:  i,  j,  p,  q  are  positive  integers  ranging  from  1  to  M=N,  h  of  Eqns  (B1)-(B5)  is  the 
parameter  defined  in  Eqn  (1),  shell  planform  rxs. 

From  Eqn  (23a): 


\  s 

fe  +  —iTdy-ij  j(l  -  cos  m  cos  pit)  =  0 
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From  Eqn  (23c): 


_]2£«  -^kjVd,  -S^K^jVd, 

rs  rs  s  s 

-^md,  -— 

S  S  S  /*3 

-^KHVe, 

J!fk,iiVe, -^k,HjVe,)+ 

rs  rs 

+^k,HVc,  +^K^jVc,  +^jVc,  +^kjVc,  +^k,'fK.% 

T  S  ^  ^  ^ 

■  '^"‘‘\vf%  \ + ^k/%  %  +^k,VK  %  + 


+^^*,1*4  +^ijVdf  +^^t,(/W,  +-^k„^ljVd^)+ 

/“  r  rs  » 

+^jKd^  +^kjKd^  +^kjKe JjicosiK  cospic  -1)+ 

+^jTidy  +^kjKd^  +^kJneJ^icosjTicosgK  -1) 


=  load 
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From  Eqn  (23d): 


.  J  S  S  ^  iJ  *3 

-^Lk,Vjx%  -^Km^d, -^k,Hin% 

T  S  VS  rS  fo  to 

-6D„t.«s -6F„K^e,  -^iVe,  -^kjVe^ 

-S^k,^jVe^)+ 

s  s  s 

(fJit-^kjVcn  -S^k,V%%  -^kMn  -^kjjVc, 

-^k^Vc,  -S^kHiVc,  -^k,HjVc,  -A,,d, -6D„k,d,-9F^k,% 

rs  rs  rs 

-B^iVd,  -^k^Vd, -i^k.Wd,  -^jVd,  -^kjVd, 

T  T  T  S  A 

5*  f^S  I A  '*3  ^*3  '*3 

-"kkk;ijn%)  + 
rs 

(/g/^^iyjc^  jnd^  +^kjndy  +^iKeg  +-^^*i7ceJ(l-cosz7ccosp7c)+ 

\  rs  s  s  r  r  j 


k/^jn^Cy  jKdy  +^kjndy  +^ii:ey  +^kjiiey  k*(l-cosijccosp7c)=0 
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From  Eqn  (23e): 


I  j  S  S  o  o  ^ 

-^k,VjK% 

1“  s  f  S  r  s  r  o  ' 

-^ijn% -‘^k.m^d^  -^kJlM%  -—k"iin% -S^k,HjK% 

rs  rs  rs  rs  rs 

-A^e, -6D„ke, -9F^k\  -^iVef-^UVe,  -i^kV^e, 

-^jVe,  -^kjVe,  -Hfk,^  iH\)A 

s  s  s 


3/^26  I  •  *2  3  ^^26  7  2  •  *2  3 

— fKij^n^c^ - ij  n 

rs  rs 


- ^Asdij  -^D^sKdy  -9F4sK^d^ 


-^kfK%  -^k,Wdy  -^f%^d^  -^kjh^d^ 
r  r  r  s  s 

s  ^  rs  rs 


{fAi^^k^ijn^Ci:  +^jndy  +^kj'Kd^  +^iTiey  +^kJiiey\l-cosjncosqn)+ 

\  rs  s  s  r  r  J 

+^j%d.  +^^kjizd.  +^iKe^  +^^kjKe\^{l-cosjiicosqn)=0 
\  rs  s  s  r  r  ) 


where  in  Eqns  (B  1)-(B5): 


2  r  i-P 


0  if  i^p 


if  j  =  Q 


0  if  j^q 
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2ri 


/3=‘ 

0 


'  2rp 

u=- 

0 


2sj 

/5=' 

0 


'  2sq 

h=' 

0 


if  i  +  p  =  odd 
if  i  +  p  =  even 


if  i  +  p  =  odd 
if  i  +  p  =  even 


if  j  +  q  =  odd 
if  j  +  q  =  even 


if  j  +  q  =  odd 


if  j  +  q  =  even 
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Nonlinear  terms 

The  additional  nonlinear  terms  below  in  Eqns  (B6)-(B8)  are  included  on  the  left  hand  side 
of  Eqns  (B1)-(B3).  The  nonlinear  Aie  and  A26  terms  are  not  shown  and  i,  j,  k,  I,  m,  n,  p,  q 
are  positive  integers  ranging  from  1  to  M=N. 

FromEqn  (23a): 


'm\(kKy  r  (  1 


j  k  I 


(i4j2  +  Agg  ) 


1  11^(1  1  1  1 


2n[xy  X2  X3  x^JlK^y^  y^  >^4 


>Y«Ell-Lf_J-+±-±+J-Wi-+±+J-+i 


rs  J  2jc  I  Xj  x.^  x^ )  2k  yyi  ^3  ^4 


jK'^(kK'\  r  f  1 


r  J  2k  {  Xi  X2  X3  X4  J  2k  {yi  3^2  >'3  ^4 


A,,  f iK^f kK\  s  (  1  .  1 


2  ^  r  A  r  ;27cI>’i  y^ 


- - ^  (1-COS/7C  COStoCOSpTC)) 

3^3  yAj 


From  Eqn  (23b): 
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FromEqn  (23c): 


^ (  A  ^  1  1  1  1  1  1 

rr  *  / ‘'vtJItJ 

2AJ^¥^Wi-+±+±^±Ur±^±^A^±l_ 

W  A  y2jci^;ci  Xj  X3  xj2it{y^  y^  y^  yj 

.,fiiYri!LW±^±.±_±w±^±_±_±lu 


rjyrJ2it\xi  X2  X3  x^J2n[y^  y^  >3  y^ 


i  j  k  I 


.  f^Yi^Li-LlJ-  1 _ I _ LlJ_(J_  J__J___i 

ijki  \  ^  Jl  A-  J  27cUi  ^  X2  X3  x^  j27tUi  ’’’jj  >^3  y^ 


2aJ!^]Me1]jl±,±^±^±U±^±,±,±]. 

Vrjy^  rs  j2Kyx,  x^  x^  xj2n{yi  y^  y^  y^) 

.  r^Yriii^r±,±.±-±ur±^±.±.±l')^ 


s  j2K{xy  x^  X3  X4j2K\y^  y^  y^  y^ 


j  *  I 


4-rAYiE'iJ_(±+±-±.±iJ.(±^.±+±+i 
2/?  J  A  5  y  2jc  I  a:i  X2  X3  X4)2n[yi  y^  ys  y4 


R  \  r  J  2;c  I  Xj  x 


'^Y-Lf  J_+ J _ I _ L1j_(±+J _ I _ L 


X4j2ii[yi  y2  >'3  3'4 


2R\  r  J\r  J2k  [x^  x 


ikf^Y”L'|Jl[±+±+±+±|J_[±+±_J__i 


^4  j  27:1^1  ^2  3^3  ^4 


A2  (iT^y  rfll  1  IWfl  1  1  1 


R  \  s  J  27C  I  X 


X4  1 231  [y^  y^  y^  y^ 
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i  j  k  I  m  n 


(B8) 


where  Xi  and  y,-  are  defined  similar  to  the  of  Appendix  A.  Instead  of  indices  k,ij  for  the 
ti,  the  Xi  are  given  by  substituting  p,i,k  respectively.  The  y,  are  defined  likewise  by 
substituting  indices  q,  j,  /  respectively.  The  g,  and  the  hi  are  defined  from  the  «,•  of 
Appendix  A.  Instead  of  indices  ij,k.l  for  the  n,,  the  g,  are  given  by  substituting  i,k,m,p 
respectively.  The  hi  are  defined  likewise  by  substituting  indices  respectively. 


APPENDIX  C 
FORTRAN  LISTING 
SAMPLE  INPUT/OUTPUT  LISTINGS 


C  21mar95,  galeikm  JOR/even  and  odd  terms 

C  GALERKIN  SOLUTION  TO  DONNELL  SHALLOW  SHELL,  LAMINATED  MATERIAL 
C  PARABOLIC  TRANSVERSE  SHEAR  residual  force 
C 

C  to  diange  to  odd  terms  only,  look  for  C  $$$$$ 

C  change  to  i=l4iterms*2-l,2 

0(i**************,*******+*************************************** 

C  INPUT: 

C  NANAL(l)  (ISO=l,  LAM=2),  NANAL(2)  (LINEAR=1,  NONLINEAR  =2) 

C  FOR  NONLINEAR:  NINCUMAX 
C  NTERMS,  count  only  odd  terms 
C  NPRINT 

C  NLOAD(SIN=0.  UNIFC«M=I,  CENTER  PT  LOAD=2),  MAG  LOAr>=QO 

C  EYJTOJIT  *OR*  E1E2.G12J«J12,G13,G23 

CFORLAM:!^^^ 

CTH1,TH2,TH3...NP 

C  R,SMD  OF  ELAT  PLATE,  INPUT  RAD=0.) 

C 

C  PANEL  RXS,  RADIUS=RAD 

(>**♦*♦****♦*♦***♦***♦**♦♦♦*♦♦*♦*♦**********♦*******♦*♦****♦**** 

c  _ 

C  ARRAYS  -STTPn  FOR  21 TERMS-LINEAR  AND  15  TERMS-NONLINEAR 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

CHARACTER*64  FNAME,GNAME 

COMMON/STF/R,SJRAD,ST(2205,2205), NTERMS  J'LOAD(2205),QO,NLOADPHO, 

X  STNL13(225,225),STNL23(225^25),STNL31(225,225). 

X  STNL32(225,225),STNL33(225,225),SnFF(2205,2205) 
COMMON/ELAS/NANAL(3)E132,G12J4U12J>T,THE(20)J4P,A(33). 

X  B(3,3)T>(3,3)E(3,3)E(3,3)JI{33), 

X  AS(22)J)S(2,2)ES(2.2),G13, 

X  G23EYJWjn'RTHE(20),GS 

COMMON/SOLVE/INDX{2205),W(2205) 
COMMON/CONVERGE/DIS(2205)4)ISPREV(2205) 
COMMON/STRS/QllTOP,Q12TOP,Q22TOP,Q16TOP,Q26TOP, 

X  QS11,QS22,Q12,Q22 

DOUBLE  PRECISION  K1,NUJ1U12JW21 

WRnE(*,923) 

READ(*,925)FNAME 

WRnE(*,924) 

READ(*,925)GNAME 

923  FORMATC  WHAT  IS  YOUR  INPUT  FILE  NAMET) 

924  FORMATC  WHAT  IS  YOUR  OUTPUT  FILE  NAME?") 

925  FORMAT  (A) 

OPEN(5JTLE=FNAME) 

OPEN(6EILE=GNAME) 

nER=0 

KTINC=1 

CALL  INPUT(NPRNTJflNC4MAX) 

CALL  ELAST(NPRNT,SPHO) 

CALL  STIFFNESS 
20CALLNLST(nER) 

CALL  SUBST(NANAL(2)4TER) 

CALLLUDCMP 

CALLLUBKSB 

CALL  CONVRGE(NCVJIERJMAX,KTINQ 
CALL  DISPLACEMENT(NTERMS  JJCVjrr4’HOR,S24ANAL(2)) 
IF(KTINCLENINC  JU4D.  NANAL(2)  .EQ.  2)GOTO  20 
END 
C  19oct94 

c 

SUBROUTINE  STIFFNESS 
C 

. . 

C  THIS  SUBROUTINE  FILLS  LINEAR  STIFFNESS  MATIUX 

C  Kd=f 

. . 

IMPLICITDOUBLE  PRECISION  (A-H,0-Z) 

COMMON/STF/R,SJL\D,ST(2205^205).NTERMSJ>LOAD(2205),QO.NLOADJ'HO, 


X  STNL13(225,225),STNL23(225^25),STNL31(225;Z25^ 

X  STNL32(225,225),STNL33(225,225),STIFF{2205^205) 
COMMON/ELAS/NANAL(3)^1^2,G12JW12JT,THE(20)»NP,A(3,3), 
X  B(3,3)J)(3,3)^(3,3)J^(3,3)JI(3^), 

X  AS(2^)J)S(2^)JS(2,2),G13, 

X  G23^YaWjrrjtTHE(20),GS 

C 

DOUBLE  PRECISION  K1,NUJ^U12^21 

PI=3.14159265359 

PI4=PI*PI*PI*PI 

Kl=-4./(3*HT**2) 

M=1 

N=1 

DO  100n=U*NTERMS*NTERMS 
DO  200  JJ=U*NTERMS*NTERMS 
200ST(njJ)=:0. 

100PLOAD(II)=0. 

D0  300LP=1J^TERMS 
DO300LQ=l,NTERMS 
DO350I=l,NTERMS 
DO350J=l,NTERMS 
C 

C  SET  UP  PARAMETERS  F  FOR  INDICES  MJLP 
C 

ITESTl=MODa+LP^) 

nEST2=MOD(J+LQ^) 

IFa£QI-P)THEN 

Fl=R/2. 

F3=0. 

F4=0, 

ELSE 

F1=0. 

IF  (ITEST1£Q.0)  THEN 
F3=0. 

F4=0. 

ELSE 

F3=(2.*I*R)/(PI*a*I-LP*LP)) 

F4=(2.*LP*R)/(PI*(LP*LP-I*D) 

ENDIF 

ENDIF 

C 

C  SET  UP  PARAMETERS  FOR  INDICES  NXQ 
C 

IF(J£Q.LQ)THEN 

F2=S/2. 

F5=0. 

F6=0. 

ELSE 

F2=0. 

IF(ITEST2£Q.0)THEN 

F5=0. 

F6=0. 

ELSE 

F5=(2.*J*S)/0PI*(J*J-LQ*LQ)) 

F6=(1*LQ*S)/(PI*(LQ*LQ-J*J)) 

ENDF 

ENDIF 

C 

C  FILL  IN  STIFFNESS,  ST,  FOR  BC-1 
C 

CEQNl 

C 

A1  =:F3*F6*2.*A(1,3)*I*J*PI*PI/(R*S)+ 
XFl*F2*(A(l,l)*I*I*PI*Py(R*R)  + 

XA(33)*J*J*PI*PI/(S*S))- 

XF6*A(1.3)*J*P]yS*(l.-DCOSa*PD*DCOS(LP*PD) 

C 

B1  =  F3*F6*(A(1,3)*I*I*PI*PVCR*R)  + 

XA(2,3)*J*J*PI*PV(S*S))+ 

XFl*F2*(A(l,2)+A(3,3))*I*J*PI*Pm^S)- 


XF6*A(l,3)*I*PI/R*(l.-DCOSa*PI)*DCOS0J»PI)) 

C 

Cl  =F3*F6*A(23)*J*PI*PH0/S  + 
XF1*F2*A(U)*I*PI*PH0/R 
C 

CEQN2 

C 

A2 = F1*F2*(A(U)+A(3,3))*I*J*PI*PI/(R*S)+ 
XF4*F5*(A(1.3)*I*I*PI*PI/(R*R)+ 

XA(23)*J*J*PI*P1/(S*S))- 

XF4*A(2.3)*J*PI/S*(l.-DOOS(J*PI)*DCOS(LQ*P0) 

C 

B2  =Fl*F2*(A(3^)*I*I*PI*PI/ai*R)+ 
XA(2^)*J*J*PI*PI/(S*S)}+ 

XF4*F5*2*A(23)*I*J*PI*PI/(R*S)- 

XF4*A(2,3)*I*PI/R*(1.-DC0S(J*PI)*D00S(LQ*PD) 

c 

C2  :=i!l*F2*A(2;£)*J*PI*PHO/S  + 

XF4*F5»A(2^)*I*PI*PHO/R 

C 

CEQN3 

C 

A3  =F4*F6*A(2^)*J*PI*PHO/S  - 
XF1*F2*A(U)*I*PI*PH0/R 
C 

B3  =  F4*F6*  A(23)*I*PI*PHO/R  - 
XFl*F2*AGa)*J*Pl*PHO/S 
C 

C3  =F4*F6*(4  *H(2^)*K1*K1*PI4*I*J**3/(R*S**3)+ 
X4*H(13)*K1*K1*I**3*J*PI4/(R**3*S)+ 
X2.*AS(U)*I*J*PI*PI/(R*S)+ 
X12.*DS(i;2)*Kl*I*J*PI*P]/(R*S)+ 
X18*FS(U)*K1*K1*I*J*PI*PI/(R*S))+ 
XFl*F2*(-A(2^)*PHO*PHO  -  (H(1,1)*K1*K1*I**4*PI4/R**4  + 
XAS(2^)*I*I»PI*PV(R»R>+ 

X6  *DS(2^)*K1  •I*I»PI*PI/(R*R>+ 

X9*FS(22)*Kl*Kl*I*I*PI*PI/ai*R)+ 

XH(2^)*K1*K1*J**4*PI4/(S**4>+ 

XAS(l,l)*J*l*PI*Pl/(S*S>f 

X6*DS(1,I)*K1»J»J*PI*PI/(S»S>+ 

X9.*FS(l,l)*Kl*Kl*J*J*Pl‘Pl/(S*S>t 

XZ*(H(U)+Z*H(33))*K1*K1*I*I*J*J*PI4/(R*R*S*S)))- 

XF6*{Kl*LP*PiyR)*Z*H(13)*Kl*I*J/(R*S)»PI*PI* 

X(DCOSa*Pr)*IXX)S(LP*PI)-l.)- 

XF4*(K1*LQ»P]/S)*2.*H(23)*K1<'I*J/0R*S)*PI*PI* 

X(DCOS(J*PI)*DCOS(LQ»PD-1.) 

C 

D3=F4*F6*(F(2,3)*K1*J**3*PI*PI*P1/(S**3>4- 

XHC23)*K1*K1*J**3*PI*PI*PI/{S**3)+ 

XAS(l»2)*J*PI/S+6*DS(i;2)*Kl*J*PI/S+ 

X9*FS(U)*Kl*Kl*J*PiyS+ 

X3*F(1,3)*K1*I*I*J»PI*PI*PI/(R*R*S)+ 

X3*H(1,3)*K1*K1*I*I*J*PI*PI*PI/(R*R*S))- 

XF1*F2*(F(1,1)*K1»I**3*PI*P1*PI/(R**3)+ 

XH(1.1)*K1*K1*I**3*PI*PI*P]/(R**3)+ 

XASa;2)*I*PI/R+6.*DS(2^)*Kl*I*PI/R+ 

X9.*FS(2^)*K1*K1*I*PI/R+ 

XF(U)*Kl*I*J*I*PI*PI*PI/ai*S*S)+ 

X2.*F(33)*Kl*I*J*J*PI*PI*PI/ai*S*S>f 

XH(l^)*Kl*Kl*I*J*J*PI*PI*PI/ffl.*S*S)+ 

X2*H(3,3)*K1*K1*I*J*J*PI»PI*PI/(R*S*S))- 

XF6*(K1*LP*PI/R)*(F(1.3)*J/S*PI+H(1,3)*K1*J*PI/S)* 

X(DCOSa*PI)*DOOS(LP*PB-l.)- 

XF4*(Kl*LQ*PFS)*(F(2,3)*jyS*PI+H(2,3)*Kl*J*PiyS)* 

X(D00S(J*PD*E)C0S(LQ*PD-1.) 

C 

E3  =  F4*F6*ffJ(U)*Kl*I**3*PI*PI*PI/5l**3> 

XH(U)*K1*K1*I**3*PI*PI*PI/(R**3)+ 

XAS(U)*I*PI/R+6.*DS(U)*K1*I*PI/R+ 

X9*FS(U)*Kl*Kl*I»PIyR+ 

X3*F(23)*Kl*I*J*J*PI*PI*PI/ai*S*S)+ 


X3.*H(2.3)*K1*K1*I*J*J*PI*PI*PI/(R*S*S))- 

XF1*F2*(F(2^)*K1*J*J*J*PI*PI*H/(S**3)+ 

XH(2^)*K1*K1*J**3*PI**3/S**3+ 

XAS(l,l)*J*PiyS+6*DS(l.l)*Kl*J*PI/S+ 

X9*FS(1,1)*K1*K1*J*PI/S+ 

XF(U)*Kl*I*I»I*PI**3/0i*R*S)+ 

X2.*F(3,3)»K1*I*I*J*PI**3/(R*R*S)+ 

XH(U)*Kl*Kl*I*I*J*PI**3/(R*R*S>f 

XZ*H(3,3)*K1*K1*I*I*J*PI**3/(R*R*S))- 

XF6*(Kl*LP*PI/R)*a’(U)*I/R*PI+H(l,3)*Kl*I*PI/R)* 

X(DCOSa*PD*DOOS(LP*PD-l.)- 

XF4*(Kl*LQ*PI/S)*ff(2,3)*I/R*PI+H(2^)*Kl*I»PI/R)* 

X(DC0S(J*PI)*DC0S(LQ*PD-1.) 

C 

CEQN4 

C 

C4 =F3*F6*ff;(23)*Kl*J»*3*PI*PI*PI/(S**3)  + 
XH(2,3)*K1*K1*J**3*PI*PI*PI/P**3)+ 
XAS(U)*I*PiyS+6.*DS(U)*Kl*J*PI/S  + 
X9.*FS(U)*Kl*Kl*J*PiyS+ 
X3*F(1,3)*K1*I*I*J*PI*PI*PI/(R*R*S>+ 
X3*H(l,3)*Kl*Kl*I*I*J*PI*PI‘P]/a^*R*S)>+ 
XF1*F2*(F(1,1)*K1*I**3*PI*PI*P1/31**3)+ 
XH(l,l)*Kl*Kl*I**3*PI*PI*PI/(R**3>f 
XAS{2^)*I*PI/R+6  •DSe;i)*Kl*I*PI/R  + 
X9.*FS(2^)*K1*K1*I*PI/R+ 
XF(i;2)»Kl*I*J*J*PI*PI*PI/(R*S*S)+ 
X2.*F(33)*K1*I*J*J*PI*PI*PI/(R*S*S)+ 
XH(U)»Kl*Kl*I*J*J»PI*PI*PI/fll*S*S)+ 
X2.*H(3,3)*K1‘K1»I*J*J»PI*PI*PI/(R*S*S))- 
XF6*2.*F(l,3)*I*J»Kl*PI*PI/ffl*S)* 
X(l.-DCOSa*PD*DCOS(U»*PI))- 
XF6*K1*2.»H(1,3)*K1*I*J*PI*PI/(R*S)* 
X(l.-DCOSa*PD*DCOS(LP*PI)) 

C 

D4  =  F3*F6*(2  •D(l,3)*I*J*PI♦PI/(R*S)+ 
X4♦F(l,3)*Kl•I*J♦PI*Py(R*S)+ 
Xi*H(U)*Kl*Kl*I*J*Pl*PI/(R*S))+ 
XFl*F2*(AS(2^>+6*Kl*DS{2^>+9*FSe^)*Kl*Kl+ 
XD(1,1)*I*I*PI*P1/(R*R)+ 
X2.*F(1,1)*K1*I*I*PI*PI/(R*R)+ 
XH(1,1)*K1*K1*I*I*PI*PI/(R*R)+ 
XD(3,3)*J*J*PI*PI/(S*S)+ 

X2.*F(3,3)*K1  *I*J*PI*PI/(S*S>t 

XH(3.3)*K1*K1*J*J*PPPI/(S*S))- 

XF6*(D(U)*J*PI/S+F(1.3)*K1*J*PI/S)* 

X(l.-DCOSa*P0*DCOS(LP*PD)- 

XF6*K1*(F(1 ,3)*  J*PI/S+H{1  ^)*K1*J*PI/S)* 

X(l.-DCOSa*P0*DCOS(LP*PI)) 

C 

E4  =F3*F6*(AS(U)  +6  *DS(U)*K1+ 

X9*FS(U)*K1*K1+D(1,3)*I*I*PI*PI/(R*R)+ 

X2.*F(l,3)*Kl*I*I*PI*PI/ai*R)+ 

XH(1,3)*K1*K1*I*I*PI*PI/(R*R)+ 

XD(2,3)*I*J*PI*P1/(S*S)+ 

X2.*F@,3)*K1»J*J*PI*PI/(S*S> 

XH(2.3)*K1*K1*J*J*P1*P1/(S*S))+ 

XF1*F2»(D(U)*I*J*PI*P1/(R*S)+ 

XD(3,3)*I*J*PI*PI/(R*S)+ 

XZ*F(12)*Kl*I*J*PI*PI/ai*S)+ 

X2.*F(33)*K1*I*J*PI*PE(R*S)+ 

XH(U)*K1*K1*I*J*PI*PI/(R*S)+ 

XH(3,3)*K1*K1*I*J*PI*PI/(R*S))- 

XF6*(P(U)*I*PI/R+F(1.3)*K1*I*PI/R)* 

X(l.-r)COSa*PI)*DCOS(LP*PI))- 

XF6*Kl*(F(l,3)*I*PI/R+H(l,3)*Kin*PI/R)* 

X(l.-DCOSa*PD*E>COS(LP*PI)) 

C 

CEQN5 

C 

C5=Fl*F2*aJ(2;2)*Kl*J*J*J*PI*PI*PI/(S**3><- 
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XH(2^)*K1*K1*J**3*PI*PI*PI/(S**3)+ 

XAS(1,1)*J*PI/S+6.*DS(1.1)*K1*J*PI/S+ 

X9.*FS(1,1)*K1*K1*J*PI/S+ 

XF(U)*K1*I*I*J*PI*PI*PI/(R*R*S}+ 

X2.*F(3,3)*K1*I*I*J*PI*PI*P1/01*R*S)+ 

XH(ia*Kl*Kl*I*I*J*PI*PI*PI/(R*R*S)+ 

X2.*H(3,3)*K1*K1*I*I*J*PI*PI*PI/(R*R*S)>+ 

XF4*F5*(F(13)*Kl*I*I*I»PI*PI*PI/ai**3)+ 

XH(1,3)*K1*K1»I*I*I*PI*PI*P1/(R**3H 

XAS(U)*I*P]/R+6*DS(U)*K1*I*P1/R+ 

X9*FS(U)*K1*K1*I*PI/R+ 

X3*F(23)*Kl*I*J*J*PI*PI*PI/ai*S*S)+ 

X3.*H(2.3)*K1*K1*I*J*J*PI*PI*PI/(R*S*S))- 

XR*2.*F(2,3)*Kin*J*PI*PI/(R*S)* 

X(l.-DCOS(J*PD*DCOS(LQ*PI))- 

XF4*K1*Z*H(2;3)*K1*I*J*PI*PI/(R*S)* 

X(1.-DCX>S(J*PD*DC0S(LQ*PI)) 

D5  =F1*F2*(D(U)*I*J*PI*PI/(R*S)+ 
XD(3,3)*I*J*PI*PI/(R*S)+ 
X2.*F(U)*K1*I*J*PI*PI/(R*S)+ 
X2.*F(33)*K1*I*J*PI*P1/(R*S)+ 
XH(U)*Kl*Kl*l*J*PI*PI/(R*S>t- 
XH(3,3)*K1*K1*I*J*PI*PI/(R*S))+ 
XF4*F5*(AS(U)<-6*Kl*DS(U)+9.*Kl*Kl*FS(i;i)+ 

XD(1,3)*I*I*PI*P]/(R*R)+ 

X2.*F(U)*K1*I*I*PI*PI/(R*R)+ 

XH{1,3)*K1*K1*I*I*PI*P]/(R*R)+ 

XD(2.3)*J*J*PI*PI/(S*S)+ 

X2.*F(23)*K1*J*J*PI*PI/(S*S>+ 

XH(2.3)*K1*K1*J*J*PI*P1/(S*S))- 

XF4*(D(2.3)*J»PI/S+F(2.3)*K1»J*PI/S)* 

X(1.-IXX)S(J*PI)*IXX)S(LQ*PI))- 

XF4»K1*(F(2,3)*J*PI/S+H(23)*K1*J*PI/S)* 

X(l.-DCOS(J*PI)*DCOS(LQ*PI)) 


C 

E5  =  F1*F2*(AS(1,1)+6.*DS(1,1)*K1+ 

X9*FS(1,1)*K1*K1+ 

XD(3,3)*I*I*PI*PI/(R*R)+ 

X2.*F(33)*Kl*I*I*PI*PI/(R*R>i- 

XH(3(3)*K1*K1‘I*I*PI*PI/(R*R> 

Xr)(2^)*I*J*PI*PI/(S*S)+ 

X2.*F(2^)*K1  •J*J*PI*P1/(S*S>4- 

XH(2^)*K1*K1*J*J*PI*PI/(S*S)>+ 

XF4*F5*(4*F(2,3)*Kl*I*J*PI*PF(R*S>t 

XZ*H(2,3)*K1*K1*I*J*PI*P1/(R*S)+ 

XZ*D(2.3)*I*J*PI*PI/(R*S))- 

XF4*(D(2,3)*I*P1/R+F(2,3)*K1*I*PI/R)» 

X(l.-DCOS(J*PI)*IXX)S(LQ»PI))- 

XF4*K1*(F(2,3)*I*PI/R+H(2,3)*K1*I*PI/R)* 

X(l.-DOOS(J*PI)*DCOSa.Q*PD) 

C 

ST(MJ0  =  -A1 

ST(M^+NTERMS*NTERMS)  =  -B1 
ST04J4+2*NTERMS*NTERMS)  =  -Cl 
C 

ST(M+NTERMS*NTERMS  JO  =  -A2 

ST(M+NTERMS*NTERMS  Jl+NTERMS*NTERMS)  =  -B2 
ST(M+NTERMS*NTERMS  JI+2*NTERMS*NTERMS)  =  -C2 

C 

ST(M+2*NTERMS»NTERMS,N)=  -A3 
ST(M+2♦^rmRMS*^^reRMS,N+NTERMS♦NTERMS)  =  -B3 
ST(M+2*NTERMS*NTERMS  J4+2*NTERMS*NTERMS) = -C3 
ST(M+2*^^reRMS*NTERMS  Jl+3*NTERMS*NTERMS)=  -D3 
ST(M+2*mERMS*NTERMSJ4+4*NTERMS*NTERMS)=  -E3 
C 

ST(M+3*NTER1^*NTERMS  Jf+2*NTERMS*NTERMS)=  -C4 
ST(M+3*jriHlMS*NTERMS  Jf+3*NTER2vlS*NTERMS)=  -D4 
ST(M+3*NTERMS*NTERMSJ4+4*NTERMS*NTERMS)=  -E4 
C 


ST(M^*mERMS*KYER}^;S+2*mERMS*mESiMS)=z  -05 
ST(M+4*NrmRMS*NTERMSJ^+3*hnrERMS*NTER^  -D5 
ST(M-}4*mERMS*NTERMS  J^44*NTERMS*N^^  -ES 

c 

C  LOAD  VECTOR  FOR  UNIFORM  PRESSURE  LOADING 
C 

IFa^QXP  .AND.  J£QEQ  .AND.  NLOAD.EQ.l) 
XPLOAD{M+2*NTERMS*NTERMS)=Fl*F2*  16^Q0/(PI*PI*I*D 
C 

CLOAD  VECTORFOR  CENTER  PANEL  POINT  LOAD 
C 

IFa^Q  J-P  .AND.  J£QXQ  .AND.  NL0AD.EQ.2) 
XPLOAD(M+2*NTERMS*NTERMS)=QO*DSINa*Piy2)*DSINg*Pl/2) 

N=N+1 

350  CONTINUE 
M=M+1 
N=1 

300  CONTINUE 
C 

C  LOAD  VECTOR  FOR  SINUSOIDAL  LOADING 
C 

IFa^OAD£Q.O)PLOAD(l+2*NTERMS*NTERMS)=R*S*QO/4. 

C 

RETURN 

END 

C 

C3dec94 

c 

SUBROUTINE  NLSTOTER) 

C 

C  NONLINEAR  STIFFNESS  TERMS-NO  A16,  A26  TERMS 
C 

Q(^^i^i^^^|l^^,i|i^tti^tn^*’l^*******^i^****^l^**********’l^**********>^t‘**************** 

IMPLICIT  DOUBLE  PRECISION  (A-H,aZ) 

COMMON/STF/R,SJL\D,ST(2205^205),NTERMSJ>LOAD(2205),QO.NLOAD^HO. 
X  STNL13(225,225),STNL23(225;i25),STNL31(225^25), 

X  STNL32(225,225),STNL33(225,225),STIFF(2205^205) 
COMMON/ELAS/NANAL(3)^lJE2,G12J^U12JT,THE(20)a^.A(3,3), 

X  B(3,3)J)(3,3)JE(3,3)J?(3,3),H(3,3X 

X  AS(2^)J)S(2^)JFS(2.2).G13, 

X  G23£YJSfUjrrRTHE(20),GS 

COMMON/CONVERGE/DIS(2205)J)ISPREV(2205) 

DOUBLE  PRECISION  KLNUJW12^U21 
IF(NANAL(2)JBQ.1)RETURN 
IF(ITEREQ.O)RETURN 
PI=3.14159265359 
MM=1 
NN=1 
C 

DO  10LP=i;<TERMS 

IF(NIERMS.GT.9)PRINT*JJP 

DO  10LQ=1,NTERMS 

DO30I=l,NTERMS 

DO30J=lJ^TERMS 

C1=0. 

C2=0. 

A3=0. 

B3=0. 

C3=0. 

KT1=1 

P 

DO50K=l,NTERMS 

DO50L=ia'ITERMS 

C 

C  SET  UP  PARAMETERS  X,Y  FOR  INDICES 
C 

Il=MOD(LP.a-K)^) 

I2=MODgj>+a-I0^) 

I3=MOD(LP-a+K)^) 

I4=MOD(K+I+LP^) 
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IFai-EQO)THEN 

X1=0. 

ELSE 

Xl=17(LP-a-K)) 

ENDF 

IF  (I2£Q.O)  THEN 
X2=0. 

HLSE 

x2=uaJ’+a-K)) 

ENDF 

IFa3£Q.0)THEN 

X3=0. 

ELSE 

X3=17(LP-a+K)) 

ENDIF 

IF  a4JEQ.O)  THEN 
X4=0. 

ELSE 

X4=l./(LP+a+K)) 

ENDff 

C 

J1=M0D(LQ-(J-L)^) 

J2=M0D(LQ+(I-L)^) 

J3=M0D(LQ-(J+L)^) 

j4=mod(lq+j+l;i) 

IF  (Jl£Q.O)  THEN 
Y1=0. 

ELSE 

Y1=I./(LQ-(J-L)) 

ENDIF 

IF(J2£Q.0)THEN 

Y2=0. 

ELSE 

Y2=1./(LQ+(J-L)) 

ENDIF 

IF  (J3£Q.O)  THEN 
Y3=0. 

ELSE 

Y3=1./(LQ-(J+L)) 

ENDIF 

IF  (J4£Q.O)  THEN 
Y4=0. 

ELSE 

Y4=17(LQ+J+L) 

ENDIF 

C 

CEQNl 

C 

C11-A(1,1)*I*PI/R*(K*PI/R)**2*R*S/(4*PI*PD* 

X(XI-X3-X2+X4)*(YI+Y2-Y3-Y4) 

C12=(A(1^)+A(3,3))*J*PI/S*K/R*L/S*PI*PI*R*S/(4*PI*PD* 

X(X2-X3-X1+X4)*(Y1+Y2+Y3+Y4) 

C13=-A(33)*{J*ITyS)**2*K*PI/R*R*S/(4*PI*PD* 

X(X2-X3-XI+X4)*(Y1+Y2-Y3-Y4) 

C14=.5*A(l,l)*I*PI/R*K*PI/R*S/(2*PI)*(l.-DCOSa*PI)* 

XDCOS(K*PI)*DCOSai**PD)*(Yl+Y2-Y3-Y4) 

SUMI=(C11+CI2+CI3+CI4)*DIS(KT1+2*NTERMS»NTERMS) 

C1=CI+SUMI 

C 

CEQN2 

C 

C21=^A(2^)*J*P1/S*(L*PI/S)**2*R*S/(4*PI*PD* 

X(X1+X2-X3-X4)*(-Y2-Y3+Y1+Y4) 

C22=(A(U)+A(3,3))»I*PI/R*K/R*L/S*PI*PI*R*S/(4*PI*PI)* 

X0C1+X2+X3+X4)»(-Y1-Y3+Y2+Y4) 

C 

C23=-A(3,3)*L*PI/S*a*PI/R)**2*R*S/(4*PI*PD* 

X(X1+X2-X3-X4)*(Y2-Y3-Y1+Y4) 

C24=.5*Aa^)*I*PI/S*L*PI/S*R/(2*PD*(l.-DCOS(J*PI)* 

XDOOS(L*PI)*DCOS(LQ*PI))*(Xl+X2-X3-X4) 

SUM2=(C21+C22+C23+C24)*DIS(KT1+2*NTERMS*NTERMS) 


C2=e2+SUM2 


c 

CEQN3 

C 

A31=-A(i;Z)*(L*PiyS)**2*I*PI/R*R*S/(4*PI*PI)* 

X(X1+X2.X3-X4)*(Y1+Y2-Y3-Y4) 

A34^2*A(3,3)*K/R*LyS*PI*PI*J*PI/S*R*S/(4*PI*PD* 

X(X1+X2+X3+X4)*(Y1+Y2+Y3+Y4) 

A36=-A(1,1)*(K»P]/R)*»2»I*PI/R*R'‘S/(4*PI*PI)* 

X(X1+X2-X3-X4)*(Y1+Y2-Y3-Y4) 

SUM3=(A31+A34+A36)*DIS(KT1+2»NTERMS*NTERMS) 

A3=A3+SUM3 

C 

B31=-A(U)*(K*PI/R)**2*J*PI/S*R*S/(4»PI*PD* 

X(X1+X2-X3-X4)*(Y1+Y2-Y3-Y4) 

B34=-2‘A(33)*K/R*L/S*PI*PI*I*PI/R*R*S/(4*PI*PI)* 

X0C1+X2+X3+X4)*(Y1+Y2+Y3+Y4) 

B36:^A(2^*(L*PI/S)**2*J*P]/S*R*S/(4*PI*PI)* 

X(X1+X2-X3-X4)*(Y1+Y2.Y3-Y4) 

SUM4=(B31+B34+B3Q*DIS(KT1+2*NTERMS*NTERMS) 

B3=B3+SUM4 

C 

C31=-A(2;!)*PHO/2*J*PI/S*L*PI/S*R*S/(4*PI»PD* 

X(X1+X2-X3-X4)»(Y1+Y2+Y3+Y4) 

C32=^A(U)*PHO*(K*PI/R)**2*R*S/(4*PI*PI)* 

X(X1+X2-X3-X4)*(Y1+Y2-Y3-Y4) 

C33=-A(U)*PHO/2*I»PI/R*K*P]/R*R*S/(4*PI*PD* 

X(X1+X2+X3+X4)*(Y1+Y2-Y3-Y4) 

C39=-A(2^)*PHO*(L*PiyS)**2*R»S/(4*PI*PI)* 

X(X1+X2-X3-X4)*(Y1+Y2-Y3-Y4) 

SUM5={C31+C32+C33+C39)*DIS(KT1+2*NTERMS»NTERMS) 

C 

KT2=1 

SUM=0. 

DO70M=l,NTERMS 

D070N=1JOTERMS 

C 

XX1=0. 

XX2=0. 

XX3=0. 

XX4=:0. 

XX5=0. 

XX6=0. 

XX7=0. 

c 

IF(a-K>(M-LP)£Q.O)  XXl=R/8 
IF(a-K)-(M+LP)JEQ.O)  XX4=:R/8 
IF((I+K)-(M-LP)£Q.O)  XX6=R/8 
IF(a+K)-(M+LP)£Q.O)  XX3=R/8 
IF(g.K)+(M-LP)£Q.O)  XX2=:R/8 
F(a-K)+(M+LP)£Q.O)  XX5=R/8 
IF(a+K)+(M-LP)£Q.O)  XX7=R/8 
C 

YY1=0. 

YY2=0. 

YY3=0. 

YY4=0. 

YY5=0. 

YY6±K). 

YY7=0. 

C 

IF((J-L)-(N-LQ)^Q.O)  YYl=S/8 
IF((J-L)-(N+LQ)£Q.O)  YY5=S/8 
IF((J+L).(N-LQ)£Q.O)  YY3=S/8 
IF((J+L)-(N+LQ)£Q.O)  YY7=S/8 
IF((J-L)+(N-LQ)£Q.O)  YY2=S/8 
IF((J-L)+(N+LQ).EQ.O)  YY6=S/8 
IF((J+L)+(N-LQ)£Q.O)  YY4=S/8 
C 

C34=0.5*A(2^)*J*PI/S*L*PI/S*(N*PiyS)**2* 


X(XX1+XX2+XX3-XX4-XX5-XX6-XX7)*(YY1+YY2+YY3+YY4-YY5-YY6-YY7) 

C35=(A(U)/2)*(J*PI/S)*»2*K*P1/R*M*PI/R* 

X(XX1-XX2+XX3-XX4+XX5-XX6+XX7)*(YY1+YY2-YY3-YY4-YY5-YY6+YY7) 

C36A^(A(3,3))*I*PI/R»L*PI/S*M/R*N/S*PI*PI* 

X(XX1-XX2+XX3-XX4+XX5-XX6+XX7)*(YY1-YY2-YY3+YY4-YY5+YY6+YY7) 

C36B=-(A(3,3))*K*PlyR*J*PI/S*M/R*N/S*PI*PI» 

X(XX1-XX2+XX3-XX4+XX5-XX6+XX7)*(YY1-YY2-YY3+YY4-YY5+YY6+YY7) 

C37=(A(U)/2)*a*PI/R)**2*(L*PI/S)*N*PI/S* 

X(XX1+XX2+XX3-XX4-XX5-XX6-XX7)*(YY1-YY2-YY3+YY4-YY5+YY6+YY7) 

C38=0.5*A(l.l)*I*PI/R*K*PiyR*(M*PI/R)**2* 

X(XX1+XX2-XX3-XX4-XX5+XX6+XX7)*(YY1+YY2-YY3-YY4-YY5-YY6+YY7) 

SUM=(C34+C35+C36A+C36B+C37+C38)*DIS(KT2+2*NTERMS*NTERMS)+SUM 

KT2=KT2+1 
70  CONTINUE 

SUM6=SUM*DIS(KT1+2*NTERMS*NTERMS) 

C3=C3+SUM5+SUM6 
KT1=KT1+1 
50  CONTINUE 


C 

STNL13(MM,NN)  =  C1 
STNL23(MM,NN)  =  C2 
STNL31(MMJJN)=A3 
STNL32(MM,NN)  =  B3 
STNL33(MM.NN)  =  C3 
C 

NN=NN+1 
30  CONTINUE 
MM=MM+1 
NN=1 

10  CONTINUE 
RETURN 
END 


c 

SUBROUTINE  CONVRGE(NCV  jrER4MAX.KTINQ 


IMPLICrr  DOUBLE  PRECISION  (A-ILO-Z) 

COMMON/STE/R,S,RAD,ST(2205^205),NTERMSJ‘LOAD(2205),QO,NLOADJ>HO, 
X  STNL13(225,225).STNL23(225;225),STNL31{225225). 

X  STNL32(225,225).STNL33(225,225),STIFF(2205a205) 
COMMON/ELAS/NANAL(3)^1£2,G12;4U12J>T,THE(20)J1P.A(33), 

X  B(3,3)J)(3.3)^(3,3)J^(3.3)3{3.3), 

X  AS(2;i)T)S(2^)J=S(2,2),G13, 

X  G23^Ya4U3TJRTHECM),GS 

COMMON/CONVERGE/DIS(2205)JDISPREV(2205) 

DOUBLE  PREaSION  K1,NUJ4U12NIU21 


IF(NANAL(2)£Q.l)  NCV=1 
ff(NANAL(2)£Q.l)RETURN 
IF(nER£Q.0)THEN 
nER=l 


NCV=0 

RETURN 

ELSE 

PRINT*,’  KTINC=  '.KTINC 
PRINT*,'  nER=  'TIER 
ENDIF 
PNORM1=0. 

PNORM2=0. 

PNORM3=0. 

PNORM4=0. 

NT=5*NTERMS*NTERMS 
C  ADD  DELTA  D  TOD 
D0  5I=1,NT 

5  DIS®=DISPREV(I)+DISa) 

DO  10  I=2*NTERMS*NTERMSTIT 
PNORMl=DABS(DIS(I))**2+PNORMl 
10  PNORM2=DABS(DIS(D-DISPREV(0)**2+PNORM2 
TOL=PNORM2/PNORM1*100 
DO  30  I=1,2*NTERMS*NTERMS 
PN0RM3=DABSPIS(I))**2  +  PNORM3 


30  PNORM4=:DABS(DIS(I)-DISPREVa))**2+PNORM4 
TOLl=PNORM4/PNORM3*  100. 

PRlNT*,TOL=’.TOL 

PRINT*,TOLl=’,TOLl 

IF  (TOL  .GT.  .05  .OR.  TOLl  .GT.  1.  .OR.  ITER  1£.  2)  THEN 
rrER=nER+i 
IF(ITER.GT.IMAX)THEN 
PRINT*;  MAX  ITERATIONS’ 

STOP 

ELSE 

RETURN 

ENDIF 

ELSE 

NCV=1 

WRrrE(6,849) 

849FORMAT(//) 

WRrrE(6,850)KTINCJTER 

850FORMAT(lX;iNCREMENT=’J3;  rrERATION=  ’ J2) 

WRnE(6.849) 

C 

C  BUMP  INCREMENT  COUNT,  RESET  ITER 

C  NOTE:  NCV  IS  RESET  TO  ZERO  W/N  DISPLACEMENT 

C 

KlWOKTINC+l 

nER=l 

C  INCREMENT  LOAD 
DO20I=l,NT 

20  PLOAD(l>=PLO  AD®  *KTINa(KTINC- 1) 

ENDIF 

RETURN 

END 

c 

SUBROUTINE  DISPLACEMENT(NTERMSa^CVjrrj>HOJ^S,NANAL) 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,aZ) 

COMMON/CONVERGE/DIS(2205)J)ISPREV(2205) 

COMMON/STRS/QllTOP,Q12TOP.Q22TOP,Q16TOP,Q26TOP. 

X  QS11,QS22,Q12,Q22 

DOUBLE  PRECISION  K1,NUJW12J4U21 

PI=3.14159265359 

NT=NTERMS*NTERMS 

K1=^4./(3.*HT»*2) 

U=0. 

V=0. 

W=0. 

KT=1 

C$$$$$ 

DO  10I=1,NTERMS 
DO20J=:l,NTERMS 

W=DSINa*Pl/2)*DSIN(J*Pl/2)*DIS(KT+2*NT)+W 
U=DCOSa*PD*DSIN(J*PI/2)*DIS(KT)+U 
V=DSINa*Piy2)*DCOS(J*PD*DIS(KT+NT>V 
KT=KT+1 
20  CONTINUE 
10  CONTINUE 

PRINT*;  NCV=  ’,NC  V;  W  CENTER  =  ’,W 

IF(NCV£Q.0)RETURN 

WRrrE(6,900)W 

900  FORMATCIX,^  CENTER  = '  J)20.13) 

WRrrE(6,901)U 

901  FORMAT(lX;U(R,S/2)=:  ’4D20.13) 

WRrrE(6,902)V 

902  FORMAT(lX;V(Ry2.S)=  ’  J)20.13) 

C  CALCULATE  STRESS 

EXT=0. 

EYT=0. 

Exyr=o. 

EXYB=0. 

EXB=0. 
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EYB=0. 

ES4=:0. 

ES5=0. 

EXX=0. 

EYY=0. 

T2=HT/2 

T4=HT/4 

KT3=1 

DO30I=ljraRMS 

DO40J=l,NTERMS 

SI=DSINa*PI/2)*DSINa*Piy2) 

co=Dcosa*Piy2)*Dcosg*p]y2) 

C  STRAIN  X, Y^  @  TOP/BOTTOM  CENTER  PLATE 
EXT=.I*PI/R*PIS(KT3>T2*DIS(KT3+3*NT)-T2**3*K1* 
X(DIS(KT3+3*NT)+I*P]yR*DIS(KT3+2*NT)))*SI+EX^ 
EYTKJ*PI/S*(DIS(KT3+NT)-T2*DIS(KT3+4*NT)-T2**3*K1* 
X(DIS(KT3+4*NT)+J*PI/S*DIS(KT3+2*NT)))-DIS(KT3+2*NT)*PHO)*SI+EYT 

EXYr=(J/S*DIS(KT3)+I/R*DIS(KT3+NT)-T2*(J/S*DIS(KT3+3*NT) 

X+I/R*DIS(KT3+4*NT)>T2**3*Kl*(jyS*DIS(KT3+3*NT)+I/R*DIS(KT3+4*NT) 

X+2*I/R*J/S*DIS(KT3+2*NT)))*PI*CO+EXYT 

EXYB=(J/S*DIS(KT3)+iyR*DIS(KT3+NT)+T2*(J/S*DIS(KT3+3*NT) 

X+1/R*DIS(XT3+4*NT))+T2**3*K1*(J/S*DIS(KT3+3*NT)+I/R*DIS(KT3+4*NT) 

X+2*I/R*J/S*DIS(KT3+2*NT)))*PI*COfEXYB 

Q****************************************************^**************** 

EXB=-I*PI/R*(DIS(KT3)+T2*DIS(KT3+3*NT)+T2**3*K1* 
X(DIS(KT3+3*NT)+I*P]yR*DIS(KT3+2*NT)))*SI+EXB 
EYB^J*PiyS*(DIS(KT3+NT)+T2*DIS(KT3+4*NT)+T2**3*Kl* 
X(DIS(KT3+4*NT)+J*PI/S*DIS(KT3+2*NT)))-DIS(KT3+2*NT)*PHO)*SI+EYB 
C  SHEAR  STRAINS  @  Z=0. 

ES4=(DIS(KT3+4*NT|+J*PiyS*DIS(KT3+2*NT))* 
XDSINa*PI/2)*DCOS(J*PD+ES4 
ES5=(DIS(KT3+3*NT)+I*PI/R*DIS(KT3+2*NT))* 
XDCOSa*PD''I>SIN(J*PI/2)+ES5 
C  STRAIN  X,Y@Z=:T/4 

EXX=-I*PI/R*(DIS(KT3)+T4*DIS(KT3+3*NT)+T4»*3*K1* 

X(DIS(KT3+3*NT)+I*PVR*DIS(KT3+2*NT)))*SI+EXX 

EYYK-J*PI/S*(DIS(KT3+NT)+T4*DIS(KT3+4*NT)+T4**3*K1* 

X(DIS(KT3+4*NT)+J*PI/S*DIS(KT3+2*NT))>DIS(KT3+2*NT)*PHO)*SI+EYY 

C 

KT3=KT3+1 
40  CONTINUE 
30  CONTINUE 

C  SXT  AT  (R/2,S/2,~H/2),  S  Y  AT  (R/2,S/2,H/4),  SXY  AT  (0,03/2) 

C  SS4  AT  (R/2,S,0),  SS5  AT  (R,S/2,0) 

SXT=Q1  lTOP*EXT  +  Q12TOP*EYT  +  Q16TOP*EXYT 

SXB=QllTOP*EXB  +  Q12TOP*EYB  +  Q16TOP*EXYB 

SYT=Q12TOP*EXT  +  Q22TOP*E  YT  +  Q26TOP*EXYT 

SYB=Q12TOP*EXB  +  Q22TOP*EYB  +  Q26TOP*EXYB 

SS4=QS11*ES4 

SS5=QS22*ES5 

SY=Q12*EXX+Q22*EYY 

WRITE(6,905)SXT 

905  FORMAT(lX;SIGMA-X  AT  -Z = ^J)20.13) 

WRITE(6,906)  SXB 

906  FORMAT(lX,’SIGMA.X  AT  +Z  =  ’320.13) 

WRITE(6,907)SYT 

907  FORMAT(lX;SIGMA.Y  AT  -Z  =  ’320.13) 

WRITE(6,908)SYB 

908  FORMAT(lX,’SIGMA.Y  AT  +Z  =  ’320.13) 

WRnE(6,911)  SY 

911  FORMAT(lX,’SIGMA-Y  AT  H/4  =  *320.13) 

WRnE(6,909)  SS4 

909  FORMAT(lX;SIGMA.4  AT  (R/2,S,0)  =  *320.13) 

WRrrE(6,910)  SS5 

910  FORMAT(lX;SIGMA-5  AT  (R,S/2,0)  =  *320.13) 

C 

C  RESET  CONVERGENCE  FLAG 
C 

NCV=0 


RETURN 

END 

C4NOV94 

c 

SUBROUTINE  INPUT(NPRNTJ^C  JMAX) 

C 

Q^**********************nt******j¥********^**ii^nt***n>**************i^***^f** 

IMPLICIT  DOUBLE  PRECISION  (A-H,aZ) 

COMMON/STF/R,SJIAD,ST(2205^205),NTERMSJ>LOAD(2205),QO,NLOADJ>HO, 
X  STNL13(225,225),STNL23(225^25),STNL3I(225^25), 

X  STNL32(225,225),STNL33(225,225)^TIFF(2205^205) 
COMMON/ELAS/NANAL(3)JE1^2,G12J^U12jn:,THE(20)^,A(3,3), 

X  B(3,3)J)(3,3)^(3.3)J^(3,3)J1(3,3). 

X  AS(2^)X)S(2^)J3S(2,2),G13, 

X  G23^YJWJITJITHE(20),GS 

DOUBLE  PRECISION  KLNUJW12^21 
C  READ{5,910)TnLE 
C  910  FORMAT  (20A4) 

READ(5  ♦)  NANAL(l)  J^ANAL(2) 

IF(NANAL(2)£Q.2)READ(5  *)  NINCJMAX 
READ(5  ♦)NTERMS 
READ(5*)NPRNT 
READ(5,*)NLOAD,QO 
C 

C  ECHO  INPUT 
C 

WRnE(6,896) 

896  FORMATC//) 

WRnE(6,897) 

897  FORMAT(1X,’NANAL(1)=O.UFORARB,ISO,SYMO 
WRITE(6,895) 

895  FORMAT(1X;NANAL(2)=0,UFORNL^INJEIGENO 
WRnE(6,899)NANAL(l)J^ANAL(2) 

899  F0RMAT(1X;NANAL(1):=  Ml ;  NANAL(2)=  M2) 
IF(NANAL(2).EQ.2)WRnE(6,894)NINC 
894  F0RMAT(1X,’  NUMBER  OF  LOAD  INCREMENTS^  M3) 

WRnE(6,896) 

WRnE(6,898)NTERMS 

898  FORMAT(lX,'NUMBER  OF  TERMS  IN  SERIES  (ODD/EVEN)=  ’  J2) 
WRrrE(6.896) 

WRITE(6,893)NLOAD 

893  FORMAT(1X;NLOAD(0=SINUSOIDAL,  1=UNIF0RM,2=CTRPTL0AD)  =*J2) 
WRnE(6,892)QO 

892  F0RMAT(1X;MAGNITUDE  OF  LOAD  =  ’4)20.13) 

\VRnE(6,896) 

IF  (NANAL(1)£Q.2)  GOTO  75 
C 

C  FOR  ISOTROPIC 
C 

READ(5*)EY,NU4rr 

WRnE(6,901) 

901  F0RMAT(1X.’THE  following  PROPERTIES  WERE  INPUT  ®4W,THICK)’) 
WRrrE(6,906)EY 
WRITE(6,908)NU 
WRITE(6,906)HT 
GOTO  78 
C 

C  FOR  ORTHOTROPIC,  INPUT  MATERIAL  PROPERTIES,  E1£2,G124W12 
C 

75  READ(5,*)E1,E2,G12,NU12,G13.G23 
WRITE(6,904) 

904  FORMAT(lX,’THE  FOLLOWING  PROPERTIES  WERE  INPUT) 

WRnE(6,905) 

905  F0RMAT(1X,E1,E2,G12.NU12,G13.G23') 

WRnE(6.906)El 

WRnE(6,906)E2 

WRITE(6,906)G12 

WRnE(6,908)NU12 

WRITE(6,906)G13 

WRITE(6,906)G23 
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906FORMAT(1XJ)20.13) 

908  FORMAT(1XJ)20. 13) 

REAIX5,*)NPJT 

READ(5  *)  (THEaDJI=l  W 

WRITE(6.916) 

916F0RMAT(1X) 

WRrrE(6,918) 

918  FORMATCIX.THE  FOLLOWING  LAMINATE  WAS  INPUT) 

D0  76n=l.NP 

76  WRrrE(6,920)THE(II) 

920FORMAT(1XJ)20.13) 

WRITE(6.916) 

WRrrE(6,922)PT 

922  F0RMAT(1X,’PLY  THICKNESS  =  ’  J)20.13) 

78READ(5,*)R,SJIAD 

WRrrE(6,919)R.SMD 

IF(RAD£Q.0.)THEN 

PHO=0. 

ELSE 

PHO=l./RAD 

ENDIF 

919  F0RMAT(1X,TANEL' J)20.13,'  X ' J)20.13;  RAD=  '0)20.13) 
WRnE(6,896) 

RETURN 

END 

C 

SUBROUTINE  ELAST(NPRNT,S  J>HO) 

Q 

(^♦♦♦♦♦♦♦♦♦★♦♦♦♦★♦♦♦♦♦♦♦♦★♦♦♦★♦♦♦♦♦♦♦♦******************************** 

C  THIS  SUBROUTINE  CALCULATES  THE  ELASHCrrY  MATRICES, 

C  HOVSJSES 


IMPLICIT  DOUBLE  PRECISION  (A-H,0-2) 

COMMON/ELAS/NANAL(3),E1,E2,G12JJU12JU.THE(20)JIP,A(33), 
X  B(3,3)J)(3,3)«3, 3)05(3, 3), H(3.3). 

X  AS(2,2)0)S(2,2)05S{2,2),G13, 

X  G23,EYJ4UJITOITHE(20),GS 

COMMONySTRS/QlITOP,QI2TOP,Q22TOP,Q16TOP,Q26TOP, 

X  QSI1,QS22,Q12,Q22 

DOUBLE  PRECISION  K1,NUJIU120IU2I 
DIMENSION  QBAR(3,3),QSBAR(2,2) 

C 

C  INITIALIZE  THE  ELASTICITY  MATRICES 
C 

DO10M=I,3 
DO  il  N=l,3 
A(MJI)=0. 

B(MOI)=0. 

D(M0I)=0. 

E(MOI)=0. 

F(M0I)=O. 

11  H(MJ^. 

10  CONTINUE 
DO  15  M=U 
D016N=1,2 
AS(MOO=0. 

DS(M0r)=0. 

16  FS(M,N)=0. 

15  CONTINUE 
IF(NANAL(l).NE.l)GOTO  30 
C 

C  ISOTROPIC  CASE 
C 

GS=EY/(2.*(1.+NU)) 

DENOM=l.-NU**2 

QBAR(l,l)=EY/DENOM 

QBAR(U)=NU*EY/DENOM 

QBAR(2,2)=QBAR(1,1) 

QBAR(2,1)=QBAR(1,2) 


nnn  ^nnnn 
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QBAR(3,3)=GS 

QBAR(1.3>=0. 

QBAR(3.1)=0. 

QBAR(2,3>=0. 

QBAR(3^)=0. 

QSBAR(1,1)=GS 

QSBAR(2^)=GS 

QSBAR(U)=0. 

QSBAR(2,1)=0. 

QllTOP=QBAR(l.l) 

Q12TOP=QBAR(U) 

Q22TOP=QBAR(2^) 

QS11=QSBAR(1,1) 

QS22=QSBAR(2^) 

QS12=0. 

Q22=QBAR(2^) 

DO20M=l,3 

D0  21N=:1,3 

A(NW=QBAR(NtN)*HT 

D(\W=QBAR(MJ^)*HT**3/(3*2,**2) 

F(MJ^)<)BAR(NW*HT**5/(5*2.**4) 

21  H(MJS0==QBAR(NW*HT**7/(7*2**6) 

20  CONTINUE 
D0  25M:=U 
IX)26N=U 

AS(M^=QSBAR(KW*HT 

I>S(MJ4)<^SBAR(Ng^*HT**3/(3*2.**2) 

26  FS(MJS0==QSBAR(NW*HP'*5/(5*2.**4) 

25  CONTINUE 
GOTO  29 

CALCULATE  REDUCED  STIFFNESSES 
FOR  LAMINATED  ANISOTROPIC  STRUCTURES 

JO  HT=PT*NP 
NU21=E2*NU12/E1 
DENOM=L-NU12*NU21 
Qll=El/DENOM 
Q12=NU12*E2/DENOM 
Q22=E2/DENOM 

CALCULATE  INVARIANTS 


U1=(3.*Q1 1+3  *Q22+2  *Q12+4.*G12)/8. 

U2=<Qll-Q22)/2. 

U3=(Qll+Q22.2*Q12.4.»G12)/8. 

U4=(Q1 1+Q22+6  *012-4  ♦G12)/8. 

U5KQ1  l+Q22-Z*Q12+4.*G12)/8. 

C  CALCULATE  THE  ELASTICITY  MATRICES  * 

C  ♦ 

C  REMEM  THAT  THE  Z  AXIS  POINTS  DOWN  AS  IN  JONES  ♦ 

C  HOWEVER,  THE  FIRST  PLY  IS  THE  TOP  PLY.  IE,  * 

C  THE  PLY  WITH  THE  MOST  NEGATIVE  ZI!!  * 

D045n=l,NP 

45  RTHE(II)=THE(ID*3.14159265359/180. 

DO50KK=lJ4P 

QBAR(1.1)=U1  +  U2*DCOS(l*RTHE(KK))  +  U3*DCOS(4.*RTHE(KK)) 
QB  AR(1,2)=U4  -  U3*DCOS(4  *RTHE(KK)) 

QBAR(2,2)=U1  -  U2*DCOSC  *RTHE(KK))  +  U3*DCOS(4  *RTHE(KK)) 
QBAR(1,3)=.5*U2*DSIN(2  ♦RTHE(KK))  +  U3*DSIN(4  ♦RTHE(KK)) 
QBAR(2,3)=.5*U2*DSIN(Z*RTHE(KK))  -  U3*DSIN(4  ♦RTHE(KK)) 
QBAR(3.3)=U5  -  U3*DCOS(4  *RTHE(KK)) 

QBARa,l)=QBAR(l,2) 

QBAR(3,1)=QBAR(1,3) 

QBAR(3,2)=QBAR(2.3) 

QSBAR(l,l)=G23*DCOS(RTHE(KK))**I+G13*DSIN(RTHE(KK))**2. 

QSBAR(2,2)=G13*DCOSg^THE(KK))**Z+G23*DSIN^THEaCK))**2. 

QSBAR(1.2)^(G23-G13)*DCOS(RTHE(KK))*DSIN(RTHE(KK)) 

QSBAR(2,1)=QSBAR(1,2) 


nnn 


IF(KK.EQ.1)THEN 

QllTOP=QBAR(l,l) 

Q12T0P=QBAR(U) 

Q22TOPi=QBAR(2^) 

Q16TOP=QBAR(l,3) 

Q26TOP=QBAR(2,3) 

ELSE 

ENDF 

IF(KK.EQ.NP/2+l)THEN 

QS11=QSBAR(1,1) 

QS22=QSBAR(2^) 

ELSE 

ENDF 

1F(KK£Q.3)THEN 

Q12=QBAR(U) 

Q22=QBAR{2^) 

ELSE 

ENDIF 

ZL=(KK*1.  -  NP*^*PT 
ZU=ZL-PT 
57  DO  51  M=1.3 
D0  52N=1,3 

A(KLN)=A(M,N)  +  QBAR(M  J0*PT 
D(KLN)=D{KLN)  +  QBAR(KLN)*(ZL**3-ZU**3)/3. 
F(Ma^=F(MJ0  +  QBAR(MJl)*^**5-ZU**5)/5. 
H(\W=H(lvLN)  +  QBAR(MJ0‘(ZL**7-ZU**7)/7. 

IF  (NANAL(1)£Q.2)  GOTO  52 

B(NLN)=B(M  Jl)  +  QBAR(MJl)*(ZL**2-ZU**2)/2. 

E(MJ^:«(NLN)  +  QBAR(HN)»(ZL**4-ZU**4)/4. 

52  CONTINUE 
51  CONTINUE 
DO60M=U 
DO  61  N=U 

AS(NLN)=AS(MJ^+QSBAR(MJ^*PT 

DS(MJ4)=DS(M;4)+QSBAR(MJ4)*(ZL**3-ZU**3)/3. 

61  FS(NW=FS(NLN)+QSBAR(MJ<)*(ZL**5-ZU**5)/5. 

60  CONTINUE 
50  CONTINUE 
C 

C  SET  TO  ZERO  THOSE  ENTRIES  DUE  TO  ROUNDOFF  ERROR 
C 

29  DO  85  M=l,3 
D0  86N=1,3 

IF(DABS(A(1,1))-GT.DABS(A(NW*1.D08))A(MJO=0. 

IF(DABS(B(1,1)).GT.DABS(B(MJ4)»1.D08))B(M,N>=0. 

IF(DABS(D(1,1)).GT.DABS(D(MJ^*1.D08))D(MJ4)=0. 

IF(DABS(E(1.1)).GT.DABS(E(HN)*1.D08))E(MJJ)=0. 

IF(DABSa!(l.l)).GT.DABS(F(NW*l.D08))F(M;4)=0. 

IF(DABS(H(1.1)).GT.DABS(H(\LN)*1.D08))H(MJ^)=0. 

86  CONTINUE 
85  CONTINUE 
DO90M=U 
DO  91  N=U 

1F(DABS(AS(1,1)).GT.DABS(AS(MJ4)*1.D08))AS(NW=0. 

IF(DABS(DS(1,1))-GT.DABS(DS(M^*1.D08))DS(MJ^)=0. 

IF(DABSff^S(l,l)).GT.DABS(FS(MJ4)*l.D08))FS(Ma4)=0. 

91  CONTINUE 
90  CONTINUE 

C  CALCULATE  BATDORF-STEIN  PARAMETER 
BAT1=SQRT(A(1.1)*A(2^)-A(U)**2) 
BAT2=SQRT(12*SQRT(A(1,1)*A(2^)*D(1,1)*D(2^))) 
BATDORF=S**2‘BATl/BAT2*PHO 
WRrrE(6.962)BATDORF 

962  F0RMAT(1X,’BATIX»F  PARAMETER  =  ’4320.13) 
IF(NPRNT£Q.O)RETURN 

OUTPUT  THE  MATRICES 

WRrrE(6.916) 

916  FORMAT  (130 
WRnE(6,950) 


950  FORMATClX/AatJ^*) 

DO  65 11=1,3 

65  WRITE(6,952)A(n.l),A(n,2)AaU) 

WRrrE(6.916) 

WRITE(6,954) 

954  FORMAT(lX;Ba^=’) 

D0  66n=l,3 

66  WRrrE(6,952)B(n,l)3(n,2),B(n,3) 

WRITE(6.916) 

WRnE(6,956) 

956  F0RMAT(1X;D(I,J)=’) 

DO  67 11=1,3 

67  WRrrE(6,952)Dai,l)JD(n,2),D(n,3) 

WRITE(6,916) 

WRnE(6,958) 

958FORMAT(lX,*Eaj)=’) 

D0  68n=l,3 

68  WRITE(6,952)E(n,l)3(n,2)3(n,3) 

WRITE(6,916) 

WRrrE(6,960) 

960  F0RMAT(lX,’Fa3=') 

D0  69n=l,3 

69  WRrrE(6,952)Fai,l)JF(n,2)3(n,3) 

WRnE(6,916) 

WRnE(6,964) 

%4  FORMAT(lX;Ha,J)=') 

D07in=l,3 

71  \VRnE(6,952)H(n,l)3(n,2),Hai,3) 

\VRnE(6,916) 

WRnE(6,982) 

982  F0RMAT(lX;ASa3=  *) 

DO80n=l,2 

80  WRiTE(6,953)ASai,l)^S(n,2) 

WRnE(6,916) 

WRrrE(6,984) 

984  FORMAT(lX,’DSa,J)=’) 

DO  81 11=1,2 

81  WRlTE(6,953)DS(n,l)3)S(n,2) 

WRnE(6,916) 

WRnE(6,986) 

986  FORMAT(lX,’FSaJ)=’) 

D0  82n=l,2 

82  WRITE(6,953)FSai,l)JFSaU) 

952  FORMAT(1X,3(D20.13,2X)) 

953  FORMAT(1X,2(D20.13,2X)) 

RETURN 

END 

c 

C  SUBROUTINES  FOR  SOLVING  SIMEQ,  SOLVES  STIFF*d=DIS 
C 

SUBROUTINE  LUDCMP 

IMPLICIT  DOUBLE  PRECISION  (A.H,0-Z) 

COMMON/STF/R,S3AD,ST(2205,2205),NTERMSJ>LOAD(2205),QO,NLOAD3HO. 
X  STNL13(225,225),STNL23(225,225),STNL31(225,225), 

X  STNL32(225,225),STNL33(225,225),STIFF(2205,2205) 
COMMON/SOLVE/INDXa205).VV(2205) 

C 

N=5*NTERMS*NTERMS 

TINY=1.0D-20 

DD=1. 

DO  12 1=131 
AAMAX=0. 

DOll  J=13I 

IF(DABS(STIFFa3)-GT.AAMAX)  AAMAX=DABS(STIFFaj)) 

11  CONTINUE 

IF(AAMAX£Q.O.)  PAUSE  ’SINGULAR  MATRDC 
VV(D=L/AAMAX 

12  CONTINUE 


u  u  u 


D019J=1J^ 

IX)14I=U-1 

SUM=STIFF(IJ) 

DO  13  K=U-1 

SlM=:SUM-STIFFa«*STIFF(KJ) 

13  CONTINUE 
STIFF(I,J)=SUM 

14  CONTINUE 

AAMAX=0. 

00161=14^ 

SUM=STIFF(1J) 

D015K=:1,J-1 

SUM==SUM-STIFFaJO*STIFF(iy) 

15  CONTINUE 
STIFF(I,J)=SUM 
DUM=W(I)*DABS(SUM) 
IF(DUM.GE.AAMAX)  THEN 
IMAX=I 
AAMAX=DUM 

ENDF 

16  CONTINUE 

IFg.NE.IMAX)THEN 
DO  17  K=i;s 
DUM=STIFF(IMAXJC) 
STIFF(IMAX40==STIFFajC) 
STIFF(JJO=DUM 

17  CONTINUE 
DD=-DD 

VV(IMAX)=W(J) 

ENDIF 

INDX(J)=IMAX 

IF(STIFF(JJ)£Q.O.)  STIFF(J,J)=TINY 

IF(J.NE.N)THEN 

DUM=l./STIFF(JJf) 

DO  181=1+1^ 
STIFF(I,J)=STIFF(IJ)*DUM 

18  CONTINUE 
ENDIF 

19  CONTINUE 

RETURN 

END 


SUBROUTINE  LUBKSB 

IMPLICIT  DOUBLE  PRECISION  (A-H,aZ) 

COMMON/STF/R,SJL\D,ST(2205^205),NTERMSJ'LOAD(2205),QO,NLOADJ>HO, 
X  STNL13(225,225),STNL23(225^25),STNL31(225^25), 

X  STNL32(225,225),STNL33(225,225),STIFF(2205^205) 

COMMON/SOLVE/INDX(2205),W(2205) 

COMMON/CONVERGE/DIS(2205)J)ISPREV(2205) 

N=5*NTERMS*NTERMS 

n=o 

DO  12 1=1 
LL=INDX(I) 

SUM=DIS(LL) 

DIS(LL)=DIS(D 
IFai.NE.O)THEN 
DO  11  j=nj-i 

SUM=SUM-STIFFaj)*DIS(J) 

11  CONTINUE 

ELSE  IF  (SUM.NE.O.)  THEN 

n=i 

ENDIF 

DIS(I)=SUM 

12  CONTINUE 
DO  14  I=N,1,-1 


SUM=DIS® 

IFa.LT.N)THEN 
DO  13  J=I+1 

SUM=SUM-STIFFaJO*DIS(J) 

13  CONTINUE 
ENDIF 

DIS(I)^UM/STIFFaa) 

14  CONTINUE 
RETURN 
END 

c 

SUBROUTINE  SUBST(NANALJTER) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-IiaZD 

COMMON/STF/R,SJlAD,ST(2205a205),NTERMSJ>LOAD(2205),QO,NLOADJ>HO, 
X  STNL13(225,225),STNL23(225^25),STNL31(225^25), 

X  STNL32(225,225),STNL33(225,225),STIFF(2205^205) 
COMMON/CONVERGE/DIS(2205)J)ISPREV(2205) 

NT=NTERMS*NTERMS 
IF(NANAL£Q.l  .OR.  ITER.EQ.O)GOTO  35 
D05I=U*NT 
5  DISPREV(I)=DIS® 

35  CONTINUE 
DO10I=14*NT 
DO20J=1.5*NT 
20STIFFa,J)=STa,J) 

10DIS(I)=PLOAD® 

IF(NANAL.EQ.l  .OR.  ITER.EQ.0)  RETURN 
C  FORM  UPDATED  STIFFNESS 
DO30I=l,NT 
DO40J=l.NT 

STIFFa^♦NT+J)=STIFFa;^♦NT+J)  +  STNL13a® 
STIFF(NT+U*NT+D=STIFF(NT+U*NT+J)+STNL23a® 
STIFF(2*NT+I®=STIFF(2*NT+IJ)  +  STNL3iaj) 
STIFF(2*NT+I^+D=STIFF(2*NT+I,NT+D  +  STNL32a® 

40  STIFF(2*NT+I.2*NT+J)=STIFF(2*NT+I,2*NT+^  +  STNL33a J) 

30  CONTINUE 

C  SOLVE  FOR  RESIDUAL  FORCES,  PUT  INTO  ARRAY  DIS 
DO50I=l,5*NT 
SUMMER=0. 

DO60J=l,5*NT 

60  SUMMER=^STIFFa®*DISPREV(J)+SUMMER 
50  DISffi=SUMMER+PLOAD® 

RETURN 

END 


TNPOT: 


SAMPLE  INPUT/ODTPUT  LISTINGS 
SAMPLE  MATCHES  RESULTS  OF  TABLE  5,  =  13. 


2,1 

7 

0 

1  10000. 

21e64.e6,.5e6,.25,.5e6,.2e6 

4.. 025 
0.,90.,90.,0. 

10.. 10..100. 

OUTPUT: 


NANAL(1)=0,UFOR  ARBJSO,SYM 
NANAL(2)=0,1,2FOR  NL^IN^IGEN 
NANAL(1)=2  NANAL(2)=  1 


NUMBER  OF  TERMS  IN  SERIES  (ODD  ONLY)=  7 


NLOAD(0=SINUSOIDAL,1=UN1FORM.2=CTRPTLOAD)=  1 
MAGNITUDE  OF  LOAD  =  O.lOOOOOOOOOOOOD+05 


THE  FOLLOWING  PROPERTIES  WERE  INPUT 
E1H2,G12JIU12,G13,G23 
0.2500000000000D+08 
O.lOOOOOOOOOOOOD+07 
0.5000000000000Df06 
0.2SOOOOOOOOOOOD+00 
0.5000000000000D+06 
0.2000000000000D+06 

THE  FOLLOWING  LAMINATE  WAS  INPUT 
O.OOOOOOOOOOOOOD+00 
0.9000000000000D+02 
0.9000000000000D+02 
O.OOOOOOOOOOOOOD+00 

PLY  THICKNESS  =  0.2500000000000D-01 

PANEL  0.1000000000000D+02X  O.lOOOOOOOOOOOOD+02  RAD=  O.lOOOOOOOOOOOOD+03 


BATDORF  PARAMETER  =  0.1176984347640D+02 

NCV=  IW  CENTER  =  566.09969091488 

W  CENTER  =  0.5660996909149D+03 
U(R,S/2)=  -0.6378059133087D+00 
V(R/2,S)=  0.1809982682698D+02 
SIGMA-X  AT  -Z  =  -0.7283224645778D+08 
SIGMA-X  AT  +Z  =  0.6273132129261D+08 
SIGMA-Y  AT  -Z  =  -0.3159826271727D+07 
SIGMA-Y  AT  +Z  =  0.2666074544675D+07 
SIGMA-Y  AT  -H/4  =  -0.333 1700537210D+08 
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